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Optimal  linear  smoothing  theory  is  applied  to  the  data 
from  the  Speed  of  Sound  record  attempt  of  a  three-wheeled 
rocket  car  on  17  December  1979.  A  forward-backward  estima¬ 
tion  method  is  used  which  employs  a  seven  state  forward¬ 
running  extended  Kalman  filter  and  a  Meditch-form  backward 
recursive  "fixed-interval"  smoothing  algorithm.  Data  for 
this  analysis  is  supplied  by  a  longitudinal  accelerometer 
mounted  on  the  vehicle  and  tracking  radar  measurements  of 
range,  azimuth,  and  elevation.  States  of  interest  include 
two  components  of  vehicle  position  and  velocity,  accelero¬ 
meter  time-correlated  error,  and  radar  range  and  azimuth  bias 
errors . 

Two  iterations  of  the  forward-backward  smoothing  algo¬ 
rithm  provide  excellent  covergence  of  state  estimates  and 
error  variance.  Based  on  this  analysis  a  peak  speed  esti¬ 
mate  of  1082.028  ft/sec  or  1.008  Mach  is  obtained  at  16.85 
seconds  from  the  start  of  the  high  speed  run.  After  two 
iterations  of  the  smoother  the  standard  deviation  of  the  peak 
speed  estimate  is  reduced  to  1.055  ft/sec.  We  conclude  with 
a  confidence  level  of  nearly  one,  based  on  the  assumptions 
and  modeling  techniques  employed  in  this  analysis,  that  the 
rocket  car  did,  in  fact,  exceed  the  reference  speed  of  sound 
on  17  December  1979. 
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STOCHASTIC  ESTIMATION  APPLIED  TO  THE 
LAND  SPEED  OF  SOUND  RECORD  ATTEMPT 
BY  A  ROCKET  CAR 

I .  Introduction 

Filtering  Theory 

Any  engineering  problem  inherently  involves  the  use  of 
measured  or  calculated  information.  The  engineer  bases  his 
decisions  on  a  variety  of  issues  by  using  mathematical  models 
of  the  "real  world"  to  predict  results  obtained  by  experimen¬ 
tation.  The  mathematical  models  seldom  describe  every  factor 
impacting  a  particular  issue,  but  are  simplifications  in 
order  to  describe  the  most  important  characteristics  of  the 
problem  and  maintain  tractability .  Likewise,  no  measurement 
device  can  be  considered  "perfect",  no  matter  what  accuracy 
is  claimed.  How,  then,  does  the  engineer  meet  the  ever-in¬ 
creasing  demand  for  accuracy  when  he  must  rely  on  imperfect 
models  and  measurement  devices? 

One  method  for  obtaining  better  answers  is  to  model  the 
important  characteristics  of  a  system  and  include  the  effects 
of  model  simplifications  and  measurement  imperfections.  This 
is  the  basic  idea  behind  filtering  theory,  developed  by  sev¬ 
eral  individuals,  most  notably  Kalman.  Filtering  theory  is 
concerned  with  estimating  the  "state"  or  status  of  a  system 
of  interest  at  any  time,  t^  ,  by  incorporating  the  time 
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history  of  measurements  up  through  time  t^ 
method  of  estimation,  called  the  predictor  algorithm,  com¬ 
putes  a  state  estimate  at  any  time,  t^  ,  based  on  the  time 
history  of  measurements  received  before  time,  t ^  .  The 

"Kalman  filter"  algorithm  combines  both  prediction  of  the 
state  estimate  before  time  t^  and  correction  of  this  pre¬ 
dicted  value  based  on  the  measurements  received  up  through 
time  t^  .  With  the  advent  of  high-speed  digital  computers, 
the  "Kalman  filter"  has  proven  to  be  very  suitable  in  a  num¬ 
ber  of  applications  most  frequently  in  guidance  and  control 
problems.  In  this  implementation,  state  estimates  are  gener¬ 
ated  in  an  "on-line"  manner  utilizing  the  measurement  time 
history  up  through  the  time  at  which  an  estimate  of  state 
values  is  required.  When  the  entire  time  history  of  meas¬ 
urements  over  a  particular  time  interval  of  interest  has  been 
recorded,  "off-line"  estimation  methods  can  be  implemented. 

In  this  case,  a  "smoother"  algorithm  can  be  employed  to  gen¬ 
erate  state  estimates  based  on  all  of  the  measurements,  before 
and  after  any  time  t^  .  Since  the  smoother  algorithm  has 
more  information  available  for  state  estimation,  it  is  the 
preferred  method  for  post-run  data  analysis. 

The  Kalman  filter  algorithm  requires  a  linear  model  to 
describe  the  dynamics  of  a  particular  system  and  a  linear 
relationship  between  available  measurements  and  the  states 
of  interest.  More  often  than  not,  the  system  model  or  obser¬ 
vation  relationships  are  non-linear.  This  requires  lineariza¬ 
tion  about  some  reference  values  for  the  states  as  in  the 
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"linearized  Kalman  filter",  or  re-linearization  about  the 
current  state  estimate  provided  by  the  "extended  Kalman 
filter".  The  extended  Kalman  filter  is  most  appropriate 
when  the  nominal  (reference)  state  trajectory  is  unknown  or 
when  deviations  from  the  nominal  trajectory  may  become  severe. 

Application  of  Theory 

The  purpose  of  this  report  is  to  apply  estimation  theory 
in  post-run  data  analysis  of  the  "Budweiser  Rocket  Car"  speed 
of  sound  attempt.  The  method  used  is  to  develop  an  extended 
Kalman  filter  to  describe  the  vehicle  dynamics  and  available 
measurements  of  vehicle  position  provided  by  a  tracking  radar. 
The  events  leading  up  to  this  unique  event  have  been  well 
documented  in  many  publications,  most  notably  Road  and  Track 
Magazine,  April  1980.  It  is  not  my  purpose  here  to  provide  a 
historical  account  of  the  "Budweiser  Rocket  Car"  or  to  analyze 
Air  Force  involvement  in  this  project.  However,  some  brief 
background  information  may  be  helpful. 

Description  of  Rocket  Car  Test 

In  the  summer  of  1979  a  company  known  as  Speed  of  Sound, 
Inc.  (SOS)  requested  and  received  permission  from  the  Air 
Force  to  use  Rodgers  Dry  Lake  at  Edwards  Air  Force  Base. 

Their  purpose  was  to  drive  a  rocket-powered  land  vehicle  at 
the  reference  speed  of  sound.  If  successful,  their  rocket 
car  would  be  the  first  to  attain  such  a  speed  on  land.  The 
"Budweiser  Rocket  Car"  was  not  designed  for  high  sustained 
speeds  due  to  very  limited  fuel  storage  capability  and  short 


3 


duration  thrust  augmentation  (Sidewinder  motor).  The  SOS 
plan  was  to  reach  peak  speed  very  quickly  and  then  slow  to 
a  stop.  For  these  reasons,  an  official  "land  speed  record" 
attempt  would  not  be  made.  The  Air  Force  Flight  Test  Center 
(AFFTC)  provided  a  safe  test  area  for  conducting  high  speed 
rims.  The  speed  runs  were  not  allowed  to  impact  in  any  way 
with  normal  operations  at  the  base,  and  the  government  was 
reimbursed.  At  the  request  of  SOS,  the  AFFTC  provided  com¬ 
puter  analysis  of  the  final  high  speed  run. 

The  test  track  used  for  all  the  high  speed  runs  at 
Edwards  AFB  began  at  the  northwest  corner  of  Rodgers  Dry  Lake 
(Fig.  1.1).  The  starting  position  varied  but  was  located  ap¬ 
proximately  200  feet  from  the  lake  shore.  The  course  followed 
a  straight  line  on  a  heading  of  nearly  true  south  from  remote 
camera  site  A8  to  the  intersection  of  lakebed  runways  17/35 
and  15/33.  From  here  it  made  a  12  degree  turn  to  the  right 
and  followed  runway  17  to  the  end.  This  turn  began  at  about 
the  seven  mile  point  and  had  a  radius  of  five  miles.  Use  of 
the  curve  was  required  to  take  full  advantage  of  the  length 
of  the  lake  and  would  only  be  used  if  the  normal  deceleration 
systems  on  the  vehicle  (parachute  and  brakes)  failed.  Every 
quarter  mile  along  the  course  there  were  a  pair  of  yellow 
flags  mounted  30  feet  either  side  of  centerline  to  help  guide 
the  driver.  Special  red  flags  were  used  to  signal  the  driver 
when  to  fire  the  thrust  augmentation  system.  Approximately 
two  miles  from  the  starting  position  was  a  photoelectric 
"speed  trap"  system  erected  by  the  Federation  Internationale 
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des  Motorcyclistes  (FIM),  official  speed  recorder  for  the 
event.  The  "trap"  was  actually  a  series  of  four  precisely 
measured  gates  lined  with  lights  that  would  flash  at  very 
precise  time  intervals  to  measure  average  speed  through  the 
gate.  The  traps  were  installed  and  not  moved  during  the  runs, 
which  posed  a  difficult  problem  for  the  SOS  engineers. 

The  idea  was  to  hit  the  FIM  speed  trap  at  peak  velocity 
which  would  presumably  be  just  as  the  main  engine  or  augmented 
system  ran  out  of  fuel.  The  speed  was  gradually  raised  on 
each  successive  run  by  increasing  the  amount  of  main  engine 
fuel  and  adjusting  the  timing  for  Sidewinder  ignition.  The 
point  at  which  the  vehicle  would  reach  peak  speed  for  a  given 
fuel  load/conf iguration  was  estimated  using  a  simulation  of 
vehicle  performance.  Based  on  the  data  from  this  simulation, 
the  starting  position  was  adjusted  to  reach  peak  speed  at  the 
traps.  In  most  cases  the  peak  speed  was  seldom  achieved  at 
the  traps  but  slightly  before  due  to  underestimated  vehicle 
drag  (3).  The  speed  trap  system  proved  to  be  a  very  unreli¬ 
able  measurement  of  peak  velocity  due  to  the  problems  stated 
above . 

In  addition  to  the  FIM  speed  trap  system  several  other 
instruments  were  used  to  record  speed.  A  magnetic  pickup  on 
a  rear  wheel  was  used  to  convert  wheel  rotations  into  veloc¬ 
ity.  This  device  proved  to  be  useless  above  500  miles  per 
hour  due  to  inadequate  frequency  response  of  the  device  and  a 
severe  buildup  of  dirt  over  the  run.  A  pitot  tube  installed 
on  the  nose  of  the  vehicle  measured  air  speed.  Due  to 
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compressibility  errors  near  Mach  one  and  the  unknown  influence 
of  ground  effect  this  device  was  not  considered  reliable  (3). 
Longitudinal  and  lateral  accelerations  were  measured  by  a  set 
of  accelerometers  installed  on  the  vehicle.  Accelerometer 
data  was  not  used  to  find  velocity  by  the  SOS  engineers,  only 
as  a  check  on  the  number  of  acceleration  units  (g's)  the 
driver  was  exposed  to.  Data  from  these  devices  were  recorded 
via  frequency  modulated  (FM)  telemetry.  An  Air  Force  track¬ 
ing  radar  was  used  to  track  the  car  and  provide  a  backup  of 
vehicle  performance.  This  radar  coverage  was  considered 
training  for  the  operators  and  in  no  way  impacted  on  any  mis¬ 
sion  requirements  at  Edwards. 

On  the  final  day  of  the  high  speed  runs,  the  fuel  load 
on  the  vehicle  was  increased  to  maximum  capacity  and  a  Side¬ 
winder  motor  was  installed.  The  run  was  set  for  early  morning 
to  take  advantage  of  light  winds  and  lower  temperatures. 

Radar  coverage  was  provided  by  a  tracking  radar  located  ap¬ 
proximately  4.5  miles  from  the  starting  position  on  a  hill 
overlooking  Rodgers  Dry  Lake  (Fig.  1.1).  Temperature  at  the 
speed  trap  was  recorded  by  an  FIM  official  as  20  degrees  Fahr¬ 
enheit.  Using  the  familiar  relationship  for  the  calculation 
of  the  reference  speed  of  sound,  a, 

a  =  VyRT  (1-1) 

where 

yA  ratio  of  specific  heats  for  air  =1.4 

2  2 

R&  gas  constant  =  1715  ft  /sec  -R 

TA  temperature  in  degrees  Rankine  =  479.66  R 
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we  find 


a  =  1073.536213  ft/sec  =  931.956  mph 
The  speed  of  sound  depends  primarily  on  the  temperature  of 
the  air.  The  value  at  a  given  temperature  can  also  vary  due 
to  changes  in  relative  humidity.  For  this  analysis,  we  have 
no  information  on  the  accuracy  of  the  FIM  temperature  record¬ 
ing  system,  or  relative  humidity.  Since  the  run  was  made  in 
desert  conditions  we  assume  any  changes  to  the  calculated 
speed  of  sound  due  to  relative  humidity  can  be  ignored.  We 
also  assume  that  the  recorded  temperature  of  20  degrees  is 
exact.  The  calculated  speed  of  sound  is  used  as  a  reference 
velocity  to  compare  rocket  car  performance.  A  peak  speed  of 
931.956  mph,  therefore,  was  the  goal  of  SOS  people. 

At  0726  Pacific  Standard  time  the  main  rocket  engine 
ignited,  followed  12  seconds  later  by  ignition  of  the  Side¬ 
winder.  The  vehicle  ran  out  of  fuel  about  a  fifth  of  a  mile 
prior  to  the  speed  traps  and  thus  was  already  decelerating 
as  it  passed  through  them.  The  four  traps  showed  the  vehicle 
speed  to  be  666.234  mph,  646.725  mph,  640.112  mph,  and  632.522 
mph,  respectively  (3).  Since  the  speed  trap  measurements 
were  made  after  the  vehicle  had  reached  peak  speed  they  could 
not  be  used.  The  radar  measurements  would  have  to  provide 
the  estimate  of  top  speed.  Unfortunately,  the  radar  range 
broke  lock  at  the  critical  point  during  the  run  and  followed 
a  larger  vehicle  running  parallel  to  the  test  track  and  ap¬ 
proximately  1500  feet  beyond  the  rocket  car.  After  two 
seconds  the  radar  again  picked  up  the  car.  The  Air  Force 
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radar  became  the  only  check  of  actual  vehicle  performance 
with  the  FIM  speed  trap  system  unavailable  and  unreliable 
wheel  speed  and  air  speed  indicators.  The  radar  azimuth  and 
elevation  data  were  considered  valid  for  the  following  rea¬ 
sons  (3).  The  Air  Force  radar  operator  used  a  television 
monitor  aligned  with  the  axis  of  the  radar  dish  and  manually 
adjusted  azimuth  and  elevation  tracking  rates.  Using  a  set 
of  "cross-hairs"  on  the  monitor,  the  operator  kept  the  car 
centered  on  the  television  screen.  To  attempt  to  correct 
the  erroneous  range  data,  another  vehicle  was  driven  over 
the  tracks  of  the  rocket  car.  The  same  tracking  radar  fol¬ 
lowed  this  vehicle  and  measured  range  and  azimuth.  Azimuth 
data  from  the  rocket  car  and  this  second  vehicle  were  aligned 
and  a  corrected  range  measurement  for  the  rocket  car  was 
found.  Based  on  this  corrected  range  data,  Air  Force  com¬ 
puter  analysis  showed  three  data  points  above  the  reference 
speed  of  sound,  731.96  mph.  Speed  of  Sound,  Inc.  averaged 
these  three  points  and  claimed  a  maximum  speed  of  739.66  mph 
or  1 .0106  Mach. 

Speed  of  Sound,  Inc.  elected  not  to  make  any  more  high 
speed  runs  as  the  engineers  felt  they  could  get  no  more  per¬ 
formance  from  the  vehicle  and  the  driver  had  complained  of 
stability  problems.  Thus,  the  speed  of  739.66  mph  became  the 
"official"  figure  that  was  claimed  as  the  top  speed  of  the 
"Budweiser  Rocket  Car". 
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Method  of  Analysis 

The  data  available  for  this  analysis  includes  raw  radar 
range,  azimuth,  and  elevation  sampled  at  20  points  per  sec¬ 
ond.  Also  we  have  the  data  from  the  longitudinal  accelero¬ 
meter  which  measured  specific  force  continuously.  The  te¬ 
lemetry  data  for  the  accelerometer  is  digitized  at  250  sam¬ 
ples  per  second.  Using  these  data  sources,  a  model  for  the 
dynamics  of  the  vehicle  is  developed  and  put  in  proper  form 
for  use  in  an  extended  Kalman  filter  algorithm.  A  measure¬ 
ment  model  for  the  radar  measurements  is  used  to  relate  the 
states  of  interest  to  the  available  measurements.  To  find 
the  best  estimate  of  vehicle  performance  with  the  lowest 
achievable  error,  a  "fixed-interval"  smoother  algorithm  is 
used.  The  filtering  theory  used  in  this  analysis  and  specific 
modeling  methods  are  developed  in  the  next  two  chapters. 
Chapter  IV  will  present  the  results  of  the  extended  Kalman 
filter,  while  Chapter  V  details  two  iterations  of  the  ex¬ 
tended  Kalman  filter-fixed  interval  smoother  estimation 
scheme.  Finally,  Chapter  VI  will  present  a  hypothesis  test 
of  the  peak  velocity  estimate  and  give  a  confidence  level  for 
this  estimate. 
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II .  Background  Theory 


The  problem  as  defined  in  Chapter  I  is  to  obtain  better 
estimates  of  vehicle  position  and  velocity  by  proper  modeling 
of  vehicle  dynamics  and  measurement  devices.  By  combining 
data  from  all  measurement  sources  and  including  the  effects 
of  identifiable  errors  and  noise  through  the  use  of  a  Kalman 
filter,  one  hopes  to  get  improved  state  estimates.  The 
standard  form  for  the  model  to  describe  the  dynamics  of  a 
problem  for  which  a  Kalman  filter  is  to  be  developed  is  a 
first  order  vector  differential  equation.  Generally,  a  dis¬ 
crete-time  (sampled  data)  measurement  model  is  used  to  relate 
observations  to  the  states  of  interest.  The  basic  Kalman 
filter  equations  will  be  presented  here  with  little  explana¬ 
tion.  It  is  assumed  that  the  reader  is  generally  familiar 
with  Kalman  filtering.  An  excellent  text  on  this  subject  is 
available  by  Maybeck  (4) . 


Linear  Kalman  Filter 

The  basic  equation  to  describe  system  dynamics  has  the 
following  continuous -time  form: 


x(t)  =  f[x(t),  u(t),  t]  +  G(t)w(t)  (2-1) 

where 

x(t)  -  n-state  vector 

u(t)  -  r-input  vector  (controls) 

f  -  dynamics  vector  (possibly  non-linear) 
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G(t)  -  time  dependent  coefficient  matrix  (n-by-s) 
w(t)  -  zero-mean,  white  Gaussian  noise  s-vector  of 
strength  Q.(t)  such  that 

E[w(t)w(t+r)T]  =  Q.(t)6(T)  (s-by-s)  (2-2) 

where  6(t)  is  the  Dirac  delta  function.  Available  dis¬ 
crete-time  measurements  are  modeled  by  the  following  relation 

z(ti)  =  h[x(ti),ti]  +  y(ti)  (2-3) 

where 

t.  -  discrete  measurement  time 

1 

zft^  -  ra-vector  of  discrete  measurements 
h  -  measurement  model  vector  function  (possibly 

non-linear) 

x(t^)  -  n-state  vector 

v(tt)  -  zero  mean,  white,  Gaussian  discrete-time  m-vector 
noise  process,  independent  of  system  noise  and 
of  covariance  R(t^),  i.e.  such  that 

E[v(ti)v(t:j)T]  =  R(ti)6i:j  (m-by-m)  (2-4) 

The  initial  condition  on  the  state  is  only  Known  with 
some  uncertainty,  and  is  modeled  as  a  Gaussian  random  n-vector 
assumed  independent  of  w(t)  and  v(t^)  ,  with  mean  and 

covariance: 


E[x(t  )]  =  L 


E[(x(t0)-xj)(x(t0)-x0)T]  =  ^ 


(2-5) 

(2-6) 
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The  Kalman  filter  algorithm  is  most  easily  generated 
when  the  dynamics  and  measurement  models  are  linear  relation¬ 
ships.  If  the  vectors  f  and  h  are  linear  combinations  of  the 
states,  the  dynamics  model  and  measurement  relation  become 
linear  relationships: 

x(t)  =  F ( t)x(t)  +  B(t)u(t)  +  G(t)w(t)  (2-7) 

z(ti)  =  H(t^)x(ti)  +  v(ti)  (2-8) 

where  F(t)  and  H(t^)  become  time-dependent  (or  possibly 
time  invariant)  coefficient  matrices  of  dimensions  n-by-n 
and  m-by-n,  respectively,  and  B(t)  is  an  n-by-r  matrix 
relating  control  inputs  to  the  dynamics  model. 

The  Kalman  filter  incorporates  measurement  updates 
using  the  following  relations: 

K(ti)  =  P(ti")HT(ti)[H  (ti)P(ti“)HT(ti)  +  R(ti)]-1  (2-9) 
x(t^+)  =  x(t^“)  +  K(t^)[z^-H(ti)x(t^-)]  (2-10) 

P(t.+)  =  P(ti")  -  K(ti)H(ti)P(t.-)  (2-11) 

where 

t^“  -  before  measurement  update  at  time  t^ 

t^+  -  after  measurement  update  at  time  t^ 

K(t^)  -  Kalman  filter  gain  matrix  (m-by-n) 

A 

x(t.)  -  n-state  estimate  vector 
—  i 

-  m  vector  of  measurements 
P(t^)  -  error  covariance  matrix  (n-by-n) 
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The  state  estimate  and  covariance  are  propagated  forward  to 
the  next  sample  time  from  the  initial  condition,  x( t ^  1  +  ) 
and  P(t^_^+)  at  time  t^_ ^  ,  by  integrating 

A  A 

x(t/ti_1)  =  F(t)x(t/ti_1)+B(t)u(t)  (2-12) 

P(t/ti_1)  =  F(t)P(t/ti_1)+P(t/ti_1)FT(t) 

+G(t)Q.(t)G(t)T  (2-13) 

where  t/t^_1  ,  indicates  integration  forward  from  the  pre¬ 

vious  measurement  update  time,  t^ 

Extended  Kalman  Filter 

The  case  where  either  the  vector  of  dynamics  relations, 
f,  or  the  measurement  equation  vector,  h,  is  non-linear  in 
the  states  requires  special  consideration.  The  method  most 
commonly  used  when  system  dynamics  or  measurement  non-line¬ 
arities  exist  is  the  "extended  Kalman  filter".  The  approach 
used  in  this  method  is  to  relinearize  the  dynamics  and/or 
measurement  equations  about  the  most  recent  estimate  of  the 
state,  x(ti")  ,  at  update  time,  or  x(t/t^_^)  in  the 
ensuing  sample  period.  Thus,  the  matrices  F,  H,  K,  and  P 
are  evaluated  by  knowing  the  most  recent  estimate  of  the 
nominal  (reference)  state  trajectory. 

The  system  matrix,  F(t)  ,  in  (2-13)  and  observation 
matrix,  H(t^)  in  (2-9)  and  (2-11)  become  partial  derivative 
matrices  in  the  extended  Kalman  filter: 
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F[t;x(t/ti_1)] 


3f[x(t) ,u(t)  ,t] 
3x 


x 


A 

x(t/ti_1 ) 


(2-14) 


X  =  x(ti“)  (2-15) 

In  equations  (2-14)  and  (2-15)  the  differentiation  is  done  so 
that  the  derivative  of  a  scalar  with  respect  to  a  column  vec¬ 
tor  is  a  row  vector.  The  matrices  resulting  from  this  dif¬ 
ferentiation  have  dimensions  n-by-n  and  m-by-n,  respectively. 
These  matrices  relate  small  perturbations  in  the  state  vec¬ 
tor,  x(t),  to  changes  in  the  equations  for  x(t)  and  z(t^)  . 

The  F  matrix  is  called  the  "filter  dynamics  partial  matrix" 
and  the  H  matrix  the  "measurement  sensitivity  matrix".  De¬ 
fining  the  perturbation  of  the  state,  x(t),  from  its  current 

A 

estimate,  x(t^  )  as 

6x(t)  =  x(t)-x(ti")  (2-16) 

the  perturbation  6x(t)  is  called  the  error  state  while  x(t) 
is  the  full  state. 

We  expand  equations  (2-1)  for  x(t)  and  (2-3)  for  z(t^) 
in  a  Taylor  series  about  the  current  state  estimate  in 
powers  of  6x(t) .  Since  6x(t)  is  assumed  small,  powers  of 
ox(t)  higher  than  one  are  ignored.  We  arrive  at  the  follow¬ 
ing  linearized  perturbation  equations  in  5x( t) : 

6x(t)  =  F[t;x(t/t^_1)]6x(t)+G(t)w(t)  (2-17) 


HCti;x(ti-)]  = 


aliCx(ti)  ,ti] 


3x 


15 


(2-18) 


6z.(  t  j_ )  =  H[ti;x(ti")]6x(ti)+v(ti) 

where  F  and  H  are  defined  as  in  (2-14)  and  (2-15),  respectively. 

Equations  (2-17)  and  (2-18)  are  in  the  proper  form  for 

a 

use  in  a  conventional  filter.  Thus,  an  estimate  Sx(t.+) 

—  l 

of  the  error  state  &x{t)  can  be  made  from  perturbation 
measurements,  6z(ti)  ,  using  equations  (2-9)  through  (2-11). 
The  measurement  difference  6_z(t.;)  is  called  the  residual. 

It  is  formed  by  subtracting  the  predicted  measurements 

A 

z(t^)  from  the  actual  measurements  z(t^)  : 

6z(ti)  =  zi-h[x(ti~)  , tj^]  (2-19) 

Equation  (2-10)  provides  an  updated  estimate  of  the  error 
state,  6x(t^  )  •  By  using  equation  (2-16)  we  can  obtain 

a  + 

an  updated  estimate  of  the  whole  value  state,  x(t^  ) 

x(t/)  =  ^(ti~)+^c(ti+)  (2-20) 

This  equation  places  all  of  the  available  information  into 
the  whole-value  state  estimate.  This  allows  6x(t^+)  to  be 
reset  to  zero  for  propagation  of  the  state  estimate  to  the 
next  update  time.  At  any  time,  t^,6x(t)  has  a  conditional 

A 

mean,  6x(t.)  , 'and  conditional  covariance  Pelt.)  such 

—  1  — 0  X 

that: 

E[6x(t.)|z(ti)]  =  6x(ti)  (2-21) 

and 
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ELL5x(ti)-6x(ti)][6x(ti)-6x(ti)]1 jz(ti)]  =  P.(ti)  (2-22) 

where  Z(t^)  is  defined  as  the  entire  measurement  history- 
up  through  time  t^  .  If  we  assume  a  zero  mean  initial  con¬ 
dition  on  the  error  state. 


&(tQ)  =  0 


(2-23) 


then  i_x(t)  will  be  zero  over  the  entire  interval  between 
ipdates  such  that: 


5 x ( t )  =  0  for  t .  ,  <  t  <  t . 

—  l-l  —  —  i 


(2-24) 


With  i_x(t  ")  zero,  the  error  state  update  equation  from 


( 2-10 ) 


simplifies  to 


5x(ti+)  =  ^x(t^“)+K(t^)6z(t^)  (2-25) 


^x(ti+)  =  K(ti)  6z(ti) 


(2-26) 


which  upon  substitution  into  the  full  state  update  equation 


(2-20)  produces 


xU^)  =  x(ti-)+K(ti)6z(ti)  (2-27) 

where  5_z(t^)  is  given  by  (2-19). 

Consider  the  conditional  covariance,  P^(t^)  of  the 
error  state,  &x( t )  ,  given  by  (2-21).  We  wish  to  relate 

this  conditional  covariance  to  the  conditional  covariance, 
P(t^)  of  the  whole  state,  x(t)  .  From  (2-21) 

P^(ti)  =  Ei~[6x-6x~I[6x-6x]T|  Z  (t^  )  ] 


E[[6x][6x]T|z(ti)] 


from  (2-24) 


i 


=  e[Cx-x][x-x]T | Z ( ) ]  from  (2-16) 

=  P(ti) 

Thus,  the  error  covariance  of  x(t)  is  identical  to  that  of 
6x(t)  .  Equations  (2-13)  and  (2-11)  describe  the  propaga¬ 

tion  and  update  of  the  error  state  covariance,  P(t)  ,  with 
F ( t )  and  H(t)  in  these  equations  replaced  by  equations 
(2-14)  and  (2-15),  respectively. 

The  extended  Kalman  filter  algorithm  is  summarized  here. 
The  measurement  vector  at  time  t^,z(t^)  ,  is  incorporated 

using 

K(ti)  =  P(ti')HT[ti;x(ti')]x 

[H[ti;x(ti_)]P(ti')HT[ti;x(ti')]+R(ti)]"1  (2-28) 
x(ti+)  =  x(ti")+K(ti)[zi-h[x(ti~)  , t^^]]  (2-29) 

P(ti+)  =  P(ti-)-K(ti)H[ti;x(ti‘)]P(ti_)  (2-30) 

The  estimate  is  propagated  forward  to  the  update  time,  t^  , 
from  the  previous  update  time,  fci_l  »  by  integrating 

x(t/t^_1)  =  f[x(t/t^_1 ) ,u(t) , t]  (2-31) 

and 

P( t/ti_i )  =  F[t;x(t/ti_1)]P(t/ti_1)+P(t/ti_1)FT[t;x(t/ti_1)] 

+G(t)£(t)G(t)T  (2-32) 

from  time  t^_1  to  t^  ,  using  the  initial  conditions  pro¬ 
vided  bys 
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(2-33) 


Upon  integrating  (2-31)  and  (2-32)  to  the  next  update  time, 
ti  ,  x(t.-)  and  P(t^“)  are  defined  as 


A  _  A 

x(fci  )  =  x(ti/ti_1) 


(2-35) 


P(ti")  =  P(ti/ti_1)  (2-36) 

A 

where  the  time  notation  t^/t^_1  indicates  that  x(ti/ti_1) 
has  been  integrated  to  time  t^  but  is  conditioned  on  meas¬ 
urements  through  t^_^  only.  Thus,  the  initial  condition 
for  propagation  of  the  state  estimate  and  covariance  from 
one  sample  time  to  the  next  is  constantly  being  redeclared 

A  + 

based  on  the  most  recent  state  estimate,  x(t^_1  )  and  co- 

,  .  +  . 
variance,  ) 

F ixed-Interval  Smoother 

The  standard  Kalman  filter  algorithm  can  be  called 
"forward-running"  in  the  way  current  state  estimates  are  com¬ 
puted.  Several  authors,  including  Meditch  (5)  have  shown 
that  the  state  estimate  and  covariance  can  be  calculated 
more  accurately  by  allowing  access  to  "future"  measurements. 
The  values  for  the  state  and  covariance  at  time  t^  from  a 
forward -running  Kalman  filter  are  optimally  combined  with 
estimates  of  these  quantities  from  a  "backward-running"  Kalman 
filter.  This  type  of  algorithm  was  developed  by  Fraser  (1,2) 
and  is  called  an  "optimal  smoother". 
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The  Fraser  form  of  the  optimal  "fixed  interval"  smoother 
incorporates  a  "backward  filter"  that  propagates  the  state 
estimate  and  covariance  backward  in  time  from  the  respective 

A  + 

boundary  conditions  at  the  final  time  of  interest,  x(tf  ) 
and  P(tf+)  .  At  each  time  of  interest,  t^  ,  the  state 
estimates  and  covariances  from  the  independent  forward  and 
backward  filters  are  optimally  combined  to  yield  a  "smoothed" 
state  estimate  and  updated  covariance.  Such  an  algorithm  is 
described  in  detail  by  Maybeck  (4) .  The  optimal  smoother 
has  the  advantage  of  being  able  to  "see”  the  entire  measure¬ 
ment  history  through  the  final  time,  Z(tf)  .  Calculations 
of  the  state  mean  the  covariance  conditioned  on  all  of 
Z(tf)  ,  and  not  just  the  "current”  Z(t^)  ,  by  this  method 

can  significantly  improve  state  estimation  primarily  by 
using  the  future  data. 

Meditch  (5)  has  shown  that  a  mathematically  equivalent 
algorithm  to  the  combined  forward/backward  filter  scheme  can 
be  used.  The  Meditch  form  of  the  "fixed-interval  smoother" 
algorithm  uses  the  output  of  a  forward-running  Kalman  filter 
where  the  state  estimates  and  covariances  before  and  after 

A  _  A  +  _ 

measurement  update  x(t^  )  ,  x(t^  )  ,  P(t^  )  ,  and 

P(tA)  ,  respectively,  from  the  initial  to  final  time  have 
been  stored.  Starting  from  the  boundary  conditions  at  the 
final  time,  t^  , 

x(tf/tf)  =  x(tf+)  (2-37) 

P(tf/tf)  =  P(tf+)  (2-38) 


20 


the  smoothed  estimate  is  generated  backward  in  time  using 


x(ti/tf)  =  x(ti“)+A(ti)[x(ti+1/tf )-x(ti+1-)]  (2-39) 

where  the  "smoothing  estimator  gain  matrix"  A(t^)  is  given 
by 

A(ti)  =  P(ti+)IT(ti+1,t. )P_1(ti_1*)  (2-40) 

and  l_(t^+^,t^)  is  the  state  transition  matrix  for  prop¬ 
agating  adjoint  system  quantities  backward  in  time  (4) . 

The  covariance  of  the  zero  mean  Gaussian  estimation 
error  [x(t^)-x(t^/t„)3  can  be  generated  backward  from  the 
boundary  condition  by 

P(ti/tf)  =  P(ti+)+A(ti)[P(ti+1/tf )-P(ti+1‘)]AT(ti)  (2-41) 

The  method  of  analysis  chosen  to  analyze  available  data 
from  the  rocket  car  is  to  develop  such  a  "f ixed- interval 
smoother"  algorithm  based  on  the  state  trajectory  (time 
history)  generated  from  an  extended  Kalman  filter.  The  ex¬ 
tended  Kalman  filter  will  be  shown  as  the  appropriate  choice 
due  to  non-linear  measurement  relations.  The  amount  of 
pseudo-noise  is  adjusted  to  achieve  optimum  filter  perform¬ 
ance  in  a  "tuning"  process  described  more  fully  in  Chapter 
IV.  The  next  chapter  describes  the  methods  used  to  model 
vehicle  dynamics,  errors  to  be  estimated,  and  available  meas¬ 
urements  for  implementation  in  an  extended  Kalman  filter. 
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Ill .  Model inq  Techniques 


System  Model 

The  test  track  described  in  Chapter  I  has  no  surveyed 
positions  from  which  to  reference  vehicle  position.  The  best 
information  on  the  starting  position  is  provided  by  the  radar 
which  was  set  on  the  vehicle  for  several  minutes  before  the 
start  of  the  high  speed  run.  It  is  this  position  provided  by 
the  radar  to  which  changes  in  vehicle  position  are  referenced. 

The  coordinate  system  for  the  dynamics  of  the  rocket  car 
is  chosen  as  a  Cartesian  system  fixed  at  the  starting  point 
of  the  run.  This  system  is  shown  in  Fig.  3.1  and  has  the 
x-axis  aligned  with  the  straight  portion  of  the  test  track 
(true  south)  and  the  y-axis  aligned  with  true  east.  The  time 
interval  of  interest  is  the  first  24  seconds  of  the  run,  as 
the  vehicle  achieved  its  maximum  velocity  at  approximately 
17  seconds  into  the  run.  Thus,  the  velocity  and  position 
along  the  x  and  y  axis  is  taken  with  respect  to  a  fixed  posi¬ 
tion  on  the  earth  corresponding  to  an  inexact  starting  posi¬ 
tion.  The  elevation  of  the  car  is  ignored  due  to  minimal 
change  in  vertical  displacement  (±  20  feet).  Therefore,  the 
coordinate  frame  we  are  concerned  with  becomes  planar  or 
two-dimensional.  Post-run  inspection  of  the  test  track  in¬ 
dicated  that  the  vehicle  deviated  very  little  from  track 
centerline  (3).  Therefore,  y  components  of  position  and 
velocity  are  minimal  with  the  motion  restricted  to  the  x-axis 
almost  entirely. 
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Earth-Fixed.  Coordinate  Frame 


Rocket  Car 


Fig.  3—l(a.) .  Top  View  of  Coordinate  Systems 


Fig.  3“l(b) •  Side  View  of  Coordinate  Systems 


No  information  is  available  on  engine  thrust  or  vehicle 
drag  for  use  in  a  model  of  these  forces  as  they  relate  to 
the  acceleration  and  time-dependent  mass  of  the  vehicle.  We 
have  only  the  specific  force  measurements  from  the  longi¬ 
tudinal  accelerometer.  Thus,  the  model  chosen  to  represent 
vehicle  dynamics  is  a  straight-forward  two-dimensional  kine¬ 
matic  model.  In  such  a  model,  the  components  of  acceleration 
are  the  time  derivatives  of  the  velocity  components.  The 
velocity  components,  in  turn,  are  the  time  derivatives  of 
the  position  components.  In  addition  to  the  components  of 
position  and  velocity,  we  desire  to  model  inherent  errors  in 
the  accelerometer  and  radar.  These  error  states  are  "aug¬ 
mented"  to  the  states  of  position  and  velocity  in  order  to 
estimate  their  value  and  compensate  for  their  effect. 

No  data  from  the  lateral  accelerometer  is  available,  so 
we  cannot  estimate  accelerometer  misalignment  error.  The 
tremendous  vibration  of  the  vehicle  at  high  speed  can  con¬ 
ceivably  cause  a  time-correlated  error  in  the  longitudinal 
accelerometer.  Scale  factor  (due  to  digitizing  accelerometer 
data)  and  bias  error  are  neglected  due  to  lack  of  information 
and  observability  problems  associated  with  trying  to  estimate 
more  than  one  error  term  in  a  single  accelerometer  configura¬ 
tion. 

The  radar  range  and  azimuth  measurements  are  assumed  to 
be  corrupted  by  some  unknown  bias  (constant)  errors.  Very 
little  information  is  available  on  the  accuracy  of  the  radar 
or  types  of  inherent  errors,  so  that  these  bias  terms  are 


the  only  assumed  inaccuracies  in  these  measurements.  Thus, 
the  additional  states  to  be  estimated  include  one  error 
state  for  the  accelerometer,  and  two  error  states  for  the 
radar.  Lack  of  additional  measurements  and  resulting  ob¬ 
servability  problems  precluded  the  identification  of  any 
other  errors.  It  will  be  shown  in  the  next  chapter  that  even 
these  states  are  only  weakly  observable. 

The  states  to  be  estimated  become: 

-  position  component  along  x-axis 
X2  -  position  component  along  y-axis 
x^  -  velocity  along  x-axis 
x^  -  velocity  along  y-axis 

x^  -  longitudinal  accelerometer  time-correlated  error 
x&  -  radar  range  bias  error 
x7  -  radar  azimuth  bias  error 
The  first  four  states  are  related  by  deterministic  means 
(i.e.,  velocity  is  the  first  derivative  of  position,  etc.). 

The  error  states  are  modeled  as  stochastic  processes  in  the 
following  manner. 

The  time-correlated  error  of  the  accelerometer  is  modeled 
as  a  first  order  Gauss-Markov  process,  the  output  of  a  first 
order  lag,  which  is  driven  by  white,  Gaussian  noise.  Figure 
3.2  shows  the  output  of  the  longitudinal  accelerometer  in¬ 
dicating  the  extreme  fluctuations  in  specific  force  sensed 
by  this  device.  Such  rapid  fluctuations  are  due  to  extreme 
vibrations  caused  by  rough  ground  and  engine  "pulsing"  (3). 
Conceivably,  the  accelerometer  error  state  can  also  vary 
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due  to  such  vibration.  We  desire  the  filter  to  continually 
"expect"  changes  in  this  error  and  estimate  them.  With  no 
better  information  on  the  types  of  inherent  errors  in  this 
accelerometer,  the  first-order  lag  model  is  chosen.  The 
initial  choice  for  correlation  time,  T,  is  one  second.  This 
is  a  subjective  "guess"  of  how  much  this  error  will  vary  over 
time . 

The  error  states  on  the  radar  are  not  expected  to  vary 
significantly,  if  at  all,  during  the  24  second  time  interval 
of  interest.  Consequently,  these  states  are  modeled  as  ran¬ 
dom  walk. 

Random  walk  is  defined  as  the  output  of  an  integrator 
driven  by  white  noise,  while  a  random  constant  is  the  output 
of  an  undriven  integrator.  A  random  walk  model  for  an  error 
state  indicates  to  the  filter  that  we  are  not  "absolutely” 
sure  the  value  of  this  state  never  changes.  The  values  for 
radar  bias  errors  are  allowed  to  vary  (slowly)  over  the  time 
interval  of  interest.  By  using  such  a  model  for  the  radar 
bias  error  states,  the  filter  gains  in  these  channels  remain 
non-zero  and  so  the  filter  is  able  to  detect  changes  in 
these  states.  If  a  random  constant  model  were  used,  the 
filter  would  assume  that  once  it  calculates  a  value  for  these 
states,  they  will  not  vary.  Thus,  random  walk  is  preferable 
to  the  random  constant  model  without  detailed  information  on 
the  radar  error  characteristics. 
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The  first  order  differential  vector  dynamics  model  for 
the  propagation  of  the  seven  states  of  interest  plus  addi¬ 
tive  "driving"  noise  becomes: 

x:(t)  =  x3 
x2(t)  =  x4 
x3(t)  =  ax+x5(t) 

*4(t)  =  v4(t) 

x5(t)  =  -l/Tx5(t)+v5(t) 

x6  ( t )  =  v6  ( t ) 

x?(t)  =  w?(t) 

Here  w(t)  is  a  white  Gaussian  vector  noise  process  of 
strength  Q.(t)  over  the  time  interval  [tQ,tf]  .  Off-line 
"tuning"  of  the  system  noise  matrix,  Q.(t)  ,  can  be  used  to 

match  the  available  data  as  closely  as  desired.  This  was  not 
done  in  this  analysis.  Instead,  the  noise  matrix,  Q.(t)  is 
adjusted  "on-line"  in  a  performance  analysis  in  order  to 
achieve  lowest  possible  variance  in  the  state  estimates.  The 
results  of  this  tuning  process  are  detailed  in  the  next 
chapter . 

In  the  system  model  presented  above  the  longitudinal  ac¬ 
celerometer  output  ax  is  used  to  "drive"  the  propagation  of 
the  x  component  of  velocity.  The  output  of  the  accelerometer 
is  corrupted  by  a  time  correlated  error  which  is  expected  to 
vary  frequently  during  the  run.  With  no  lateral  accelero¬ 
meter  data  available,  we  assume  the  y  velocity  component, 
x4  ,  is  well  modeled  as  a  constant  with  zero  value.  A 
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small  amount  of  "pseudo-noise"  is  added  to  this  state  to 
allow  the  filter  to  estimate  deviations  in  the  y  velocity 
from  zero.  Small  amounts  of  pseudo-noise  are  added  to  the 
radar  error  states,  and  ,  for  the  same  reason.  The 

amount  of  pseudo-noise  is  adjusted  to  achieve  optimum  filter 
performance  in  a  "tuning"  process  described  more  fully  in 
Chapter  IV . 

The  model  for  the  accelerometer  time-correlated  error 
deserves  some  further  explanation.  We  have  chosen  a  Gauss- 
Markov  process,  which  is  the  output  of  a  first  order  lag  of 
the  form: 

x5(t)  =  -lAx5  (t)+w5  (t)  (3-1) 

where  T  is  the  correlation  time. 

The  concept  of  a  time-correlated  random  process,  x,  can 
be  generalized  from  the  "second  central  moment"  or  covariance 
matrix,  P^ft)  ,  defined  as 

Exx^)  =  E[[x(t)-mx(t)][x(t)-rn<.(t)]T]  (3-2) 

where  mx(t)  is  the  mean.  This  covariance  is  an  indication 
of  the  spread  of  values  about  the  mean  at  time  t.  General¬ 
izing  (3-2)  we  can  obtain  additional  information  about  how 
fast  x(t)  sample  values  can  change  in  time.  The  covariance 
"kernel",  £Xx^tl’t2^  »  defined  for  all  t^,  t2  in  some 
time  interval  T  as: 

Pxx(ti»t2)  =  E[[x(t1)-rsx(t1)][x(t2)-mx(t2)]T]  (3-3) 
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For  zero  mean  processes  (3-3)  becomes  the  "second  non-central 
moment"  which  generalizes  to  the  "correlation  kernel"  (4) 
defined  for  all  t^,  €T  as 

-xx(tl't2)  “  E[x(t1)x(t2)T^  (3-4) 

From  (3-3)  and  (3-4),  it  can  be  seen  that 


-xx(tl't2)  =  £xx(tl’t2)+2K(tl)^(t2)T 


(3-5) 


and  thus,  if  x(fc)  is  a  zero-mean  process,  tff(t.,t-)  = 

XX  ±  z. 

P xx(ti»t2)  •  A  scalar  time-correlated  random  process 

x(t)  of  zero-mean  has  a  correlation  kernel  function  (4)  of 


the  form: 


*xx(tl*t2) 


Pxx(tl't2)  =  C  e 


.  2  ~^tl~t2^/'r 


(3-6) 


where  T  is  the  correlation  time.  If  we  compare  this  random 
process  x(t)  to  another  zero  mean  process  y(t)  with  cor¬ 
relation  time  ten  times  that  of  x(t)  such  that 


YY 


pyy|tl't2) 


2  -[t,-t2]/10T 
o  e 


(3-7) 


we  can  say  that  there  is  higher  correlation  between  the 
values  of  y(t^)  and  y(t2)  than  between  xlt^)  and  x(t9) 
One  would  expect  a  typical  sample  of  x(t)  to  exhibit  more 
rapid  variations  in  magnitude  than  y(t) 

If  a  random  process,  x(t)  ,  is  wide-sense  stationary 
(4)  then  we  can  define  the  "power  spectral  density",  PSD,  of 
such  a  process  as  the  Fourier  transform  of  the  correlation 
function  (^(t)  where  r  is  the  time  difference  between 
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two  sample  times.  An  inverse  Fourier  transform  of  the  PSD, 
in  turn,  yields  the  "autocorrelation  function".  The  power 
spectral  density  of  a  device  to  be  modeled  can  be  plotted 
versus  frequency.  The  resulting  PSD  function  can  be  inte¬ 
grated  over  the  range  of  frequencies  to  obtain  the  mean 
squared  value  of  the  process.  This  method  of  obtaining  PSD 
and  the  related  autocorrelation  helps  describe  the  errors  in 
herent  in  a  device  to  enable  proper  modeling  of  these  errors 

If  we  had  such  PSD  information  for  the  rocket  car  ac¬ 
celerometer  we  could  make  more  realistic  modeling  decisions. 
Without  such  information  we  must  make  some  subjective  model¬ 
ing  decisions.  The  correlation  time  for  the  accelerometer 
error  model  specifies  how  much  we  expect  the  mean  squared 
value  of  this  zero  mean  process  to  vary  over  time.  Choosing 
T  as  one  second  indicates  that  the  mean  squared  value  for 
this  process  can  vary  by  approximately  63  percent  in  one 
second.  A  correlation  time  of  one  second  will  be  shown  to 
be  adequate  for  our  purposes. 

The  amount  of  driving  noise  on  the  accelerometer  error 
state,  q^  ,  can  be  obtained  from  specifying  the  desired 
steady-state  variance,  Pg(°°)  .  This  will  be  explained  in 

detail  below.  Given  a  continuous-time  system  that  can  be 
described  by  the  linear  vector  differential  equation 

x(t)  =  F(t)x(t)+G(t)w(t)  (3-8) 

where 
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E[w<t)]  =  0 

E[w(t)w(t+T)T]  =  £(t)5 (T)  (3-9) 

the  propagation  of  the  covariance  of  the  state  estimate  is 
described  bys 

P(t)  =  F(t)P(t)  +  P(t)FT(t)+G(t)Q.(t)GT(t)  (3-10) 


The  strength  of  the  driving  noise  on  the  accelerometer  error 
state,  Xj-  ,  is  found  in  the  following  manner.  Referring 
to  equation  (3-1)  we  can  solve  for  the  steady  state  value  of 
the  covariance.  In  this  case,  F  =  -1/T,  Q.  =  q,  and 
G  =  1  so  that  (3-10)  becomes: 


P5(t)  =  -2/TP5(t)+q5 

and 

P5 ( t )  =  P0^e"2(t~to)/T+q5T/2[l-e"2(t_to)/r]  ; 
For  t  =  0  we  find 


(3-11) 

t>tQ  (3-12) 


P5(to> 


o5 


(3-13a) 


and  the  steady  state  variance,  Pej(°°)  ,  becomes: 

P5  (°°)  =  q5T/2  ( 3-13b) 

By  specifying  the  steady  state  variance  we  can  solve  for  the 
amount  of  driving  noise  to  be  added  to  the  propagation  of 
the  accelerometer  error  state,  x^ 
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The  only  information  available  on  the  longitudinal  ac¬ 
celerometer  is  calibration  data  taken  before  and  after  the 
high  speed  run  of  the  rocket  car.  This  data  indicates  that 
the  accelerometer  calibration  varied  by  an  average  of  .003 
g*s  between  the  two  calibration  checks.  With  no  more  in¬ 
formation  on  the  inherent  errors  in  the  accelerometer,  we 
are  forced  to  refer  to  information  on  comparable  models.  We 
desire  a  comparable  accelerometer  which  has  p  specified  bias 
error  on  the  order  of  .003  g's.  The  model  chosen  in  this 
case  is  a  Honeywell  solid-state  low-cost  accelerometer.  The 
specified  RMS  error  for  this  model  is  listed  at  ,005g. 

Thus,  we  use  this  specified  error  as  the  steady-state  devia¬ 
tion  in  error  of  the  accelerometer  actually  used  in  this  test 
run.  Obviously,  we  need  to  vary  this  value  in  order  to  check 
filter  performance,  but  ,005g  will  serve  as  a  "first-guess". 
We  now  can  solve  for  the  driving  noise,  q5 

P5(~)  =  ( .005g) 2  =  q5T/2  (3-13c) 

and 

q5(t)  =  2 ( . 005g) 2  =  .0518  ( ft/sec2 ) 2/sec  (3-13d) 

Choosing  the  initial  variance,  Pg(0)  ,  to  be  equal  to 

the  steady-state  variance,  P^t00)  ,  results  in  a  "station¬ 

ary"  process  (4)  for  the  accelerometer  error  state, 

This  stationary  characteristic  is,  in  fact,  implemented  in 
the  system  model  for  the  rocket  car. 

The  initial  condition  of  the  state  vector  at  the  start 
of  the  run  is  known  only  with  some  uncertainty  and  is 
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modeled  as  a  Gaussian  random  variable  of  zero-mean  and 
covariance,  P  ,  such  that: 


E[x(t0)]  =5^=0 


ECCx(tn)-xo][x(to)-xo]T]  = 


£o 


(3-14) 

(3-15) 


Initially,  our  uncertainty  in  the  state  values  is  very  high. 
We  choose  the  following  values  for  P  : 


P 


— o 


10000 

10000 

ioo  o 

100 

0  .0259 

225 

. 25E-6 


(3-16) 


The  initial  guess  on  the  last  three  diagonal  elements  is 
obtained  by  using  specified  Root  Mean  Sguared  (RMS)  errors 
for  the  accelerometer,  range,  and  azimuth  respectively: 

Po5  =  ( .005g) 2  =  .0259  [ft/sec2]2 

Pq6  =  (15  ft)2  =  225  ft2 

Po7  =  (.0005  rad)2  =  .25E  -  6  rad2 

The  initial  variances  of  the  first  four  states  are  subjec¬ 
tive  values  based  on  relatively  high  uncertainty  in  state 
values.  Having  developed  the  dynamics  model  for  propagating 
the  states  of  interest,  it  is  now  necessary  to  relate  these 
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states  to  the  available  measurements  of  radar  range,  azimuth, 
and  elevation. 


Measurement  Model 

The  radar  measured  quantities  are  taken  with  respect  to 
the  coordinate  system  of  the  radar,  shown  in  Fig.  3.1.  These 
measured  values  must  be  related  to  the  vehicle  as  it  moves  in 
an  earth-fixed  coordinate  frame,  both  translated  and  rotated 
from  that  of  the  radar .  The  radar  range  measurement  can  be 
related  to  the  radar  coordinate  system  by  the  following  model: 

Range  =  (xr2+yr2+zr2 )  **  (3-17) 

where  xr>  y ^ ,  and  zr  are  distance  components  along  the 
three  axes  of  the  radar  coordinate  system.  The  radar  azimuth 
measurement  is  related  to  radar  coordinate  components  by: 

A 

Azimuth  =  arctan  (yr/xr)  (3-18) 

The  radar  measurement  of  elevation,  zr  ,  is  used  only  to 
adjust  the  value  of  range  as  the  car  proceeds  along  the  track. 
It  is  now  necessary  to  relate  radar  measured  components, 
xr*  yr»  and  zr  to  components  in  the  earth-fixed  system, 

Xj  and  *2 

The  origin  for  the  earth-fixed  coordinate  system  is 
taken  to  be  the  starting  point  of  the  run.  The  radar  is  set 
on  the  car  prior  to  the  start  of  the  run,  and  computed  values 
of  distance  Xj_  =  16104  feet,  and  yr  =  18374  feet  are 
taken  from  the  initial  range  measurement.  These  initial 
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values  are  held  constant  as  the  car  travels  down  the  track, 
in  order  to  relate  differing  range  and  azimuth  measurements 
to  changes  in  position  in  the  earth-fixed  coordinate  system. 
Labeling  the  initial  distance  along  and  y^  as  DELX 

and  DELY  respectively,  one  can  relate  the  distance  traveled 
in  the  earth-fixed  frame  to  radar  measurements  by  the  fol¬ 
lowing  translation: 


xr  =  DELX  -  y-1 

yr  =  DELY  +  x2  (3-19) 

As  the  car  moves  along  the  track,  radar  distance  xr  will 
decrease  to  a  point  where  x^.  =  0  .  At  this  point  x^  = 

DELX  and  the  azimuth  angle,  0  ,  equals  90  degrees.  Sub- 

stituting  these  translations  into  the  range  and  azimuth  equa¬ 
tions  we  find: 

Range  =  [  (DELX-x1 )  2+  (DELY+x2 )  2+zr2P  (3-20) 

and 

Azimuth  =  arctan  [ (DELY+x2)/(DELX-x1 ) ]  (3-21) 

Thus,  we  have  related  the  radar  measurements  to  the  states 
we  desire  to  estimate.  Note  that  the  180  degree  rotation  of 
the  radar  coordinate  system  to  align  with  the  earth  frame 
simplifies  the  relations.  The  measurement  model  is  summa¬ 
rized  here: 

Zi(ti)  =  [(DELX-X1(ti))2+(DELY+x2(ti))2+Zr(ti)2]*S 

+xfe  (ti)+v1  (t^  (3-22) 

36 


\ 


(3-23) 


In  these  equations  each  measurement  is  corrupted  by  an  in¬ 
herent  bias  error  state,  and  some  scalar,  zero-mean  white 
noise  to  account  for  measurement  and  model  inaccuracies . 

This  white  noise  vector  v(t^)  has  a  strength  R(t.)  during 
the  24  second  interval  of  interest. 

It  is  apparent  from  (3-22)  and  (3-23)  that  the  measure¬ 
ment  relations  are  non-linear  in  the  states.  Since  we  have 
no  reference  values  for  the  behavior  of  the  states,  espe¬ 
cially  the  error  states,  we  cannot  use  a  perturbation  model 
based  on  such  a  trajectory.  Thus,  we  choose  the  extended 
Kalman  filter  algorithm.  Applying  equations  (2-28)  through 
(2-30),  we  arrive  at  the  proper  implementation  of  the  extended 
Kalman  filter  for  the  rocket  car  analysis.  The  observation 

A  _ 

sensitivity  matrix,  H[t^;x(t^  )J  in  equations  (2-29)  and 

(2-30)  is  developed  from  (2-15).  For  convenience,  only  the 

non-zero  elements  of  this  m-by-n  matrix  are  presented  here. 

r  A  -  -i 

For  notational  ease,  HLt^;x(t^  )J  is  given  as  H  and 

a  a  _ 

x(t^  )  is  given  as  x 

H(  1 , 1 )  =  -(DELX-x1")/C(DELX-x1_)2+(DELY+x2")2+zr(ti)2]i2 
H(  1 ,2)  =  (DELY+x2“)/[(DELX-x1”)2+(DELY+x2")2+zr(ti)2]Js 
H(  1 ,6  )  =  1.0 

H(  2 , 1 )  =  (DELY+x2")/[(DELX-x1_)2+(DELY+x2-)2] 

H(  2 , 2)  =  (DELX-xi“)/[(dELX-x1")2+(DELY+x2-)2] 

H(  2,7)  =  1.0  (3-24) 
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Thus,  the  extended  Kalman  filter  incorporates  the  measure- 


ments  at  time,  t  ^ 

by  using: 

K(ti) 

=  P(ti_)HT[H  P(ti“)HT+R(ti)]'1 

(3-25) 

a  + 
x(t.  ) 

=  x(ti~)+K(ti)6z(ti) 

(3-26) 

P(ti+) 

=  P ( t i~ ) -K ( t i )H  P(ti~) 

(3-27) 

is  given  by: 

)  iZj.-Cd  (DELX-x1“)2+(DELY+x2')2+zr2]+x6"] 

_  A-  A  A 

)  z2-[arctanl (DELY+x2  )/(DELX-x1  )j+x?  ] 

(3-28) 

The  state  estimate  at  time  t.  ,  is  oroDaaated  forward  to 
the  next  sample  time  t^ 


x(t/ti_i)  = 


The  error  covariance  P(t/t^_j)  is  propagated  forward  by 
integrating  equation  (2-32)  in  which 


i-i 

y  integrating: 


x3(t/ti_i) 

x4(t/ti-l} 

x5(t/ti-1)+ax 

0 


-i/rx5(t/t,_j ) 
o 
0 


(3-29) 


where  5z(tQ 


5z(ti)  = 


6zi(ti 


522(ti 
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F[t;x(t/ti_1)] 

The  initial  value  for  the  system  noise  matrix  in  (2-32)  is 
chosen  as: 

0 

0 

.01 

0  .0518 

1 

0.1E- 

The  initial  values  for  the  diagonal  elements  of  Q.(t)  are 
chosen  to  indicate  our  relative  uncertainty  of  the  behavior 
of  the  corresponding  states  over  the  time  interval  of  inter¬ 
est.  The  initial  value  for  driving  noise  on  the  accelero¬ 
meter  error  state,  x^,  has  been  previously  calculated  (3-13d) 
The  remaining  diagonal  elements  of  Q.(t)  are  chosen  by  subjec¬ 
tively  deciding  how  much  these  states  will  vary  from  con¬ 
stant  values.  By  comparing  Q44(t),  Qgg(t)  and  Q77(t)  one 
can  see  that  we  have  little  doubt  that  the  azimuth  bias  error 
x-j,  is  a  constant.  We  are  less  certain  about  the  behavior 
of  the  y  velocity,  x^,  and  even  more  uncertain  about  the 
range  bias  error  state,  Xg .  The  initial  value  for  Q44(t) 
is  based  on  our  knowledge  of  the  rocket  car  trajectory  - 
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described  as  nearly  a  straight  line  along  the  x-axis  (3). 

We  do  not  expect  the  y  velocity  to  vary  significantly  from  a 
constant  value  of  zero.  Use  of  the  radar  in  a  "look-down" 
mode  increases  our  uncertainty  in  how  the  range  bias  state, 
Xg,  will  behave  over  the  24  second  interval.  The  azimuth 
bias  error  is  not  expected  to  vary  significantly  from  the  be¬ 
havior  of  a  constant.  The  initial  value  for  Q.(t)  presented 
in  (3-30)  is  adjusted  in  a  "tuning"  process  to  achieve  lowest 
possible  variance  values  for  the  seven  states  of  interest. 
This  performance  analysis  is  detailed  in  the  next  chapter 
along  with  the  results  of  the  seven  state  extended  Kalman 
filter  developed  here. 
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IV.  Development  and  Performance  Analysis 
of  the  Extended  Kalman  Filter 

The  previous  chapters  have  developed  an  extended  Kalman 
Filter  estimation  algorithm  for  post-run  data  analysis  of 
the  Budweiser  Rocket  Car.  This  chapter  describes  the  re¬ 
sults  obtained  by  using  the  seven  state  filter  outlined  pre¬ 
viously.  The  amount  of  computer  programming  required  by 
this  analysis  is  minimal  due  to  existing  software  available 
for  the  development  of  a  Kalman  filter. 

The  computer  software  used  in  this  analysis  is  a  Monte 
Carlo  Simulation  for  Optimal  Filter  Evaluation  (SOFE)  avail¬ 
able  at  Wright-Patterson  AFB  (6,7).  The  program  was  devel¬ 
oped  under  contract  by  the  Air  Force  Avionics  Laboratory 
(AFAL)  and  is  well  documented  by  Musick  (6).  SOFE  is  invalu¬ 
able  when  designing  Kalman  Filters.  The  normal  method  used 
is  a  Monte  Carlo  analysis  (4)  whereby  a  suboptimal  (reduced 
order)  Kalman  filter  is  evaluated  against  a  "truth  model". 

The  suboptimal  filter  is  adjusted  to  achieve  the  best  pos¬ 
sible  performance  when  compared  to  a  much  higher  order  "truth 
model".  The  idea  is  to  track  the  important  characteristics 
of  a  physical  system  adequately  using  a  simpler  model.  Such 
a  reduced  order  filter  could  then  be  implemented  in  an  opera¬ 
tional  system  where  computer  capability  may  be  limited. 

For  the  purpose  of  analyzing  the  rocket  car  data,  SOFE 
is  used  to  integrate  the  dynamics  equations  and  update  the 
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states  of  interest  at  measurement  times.  There  is  no  truth 
model  available  for  tuning  purposes,  nor  do  we  have  the 
ability  to  generate  sample  statistics.  SOFE  implements  the 
Kalman  filter  equations  for  either  the  linear  or  extended 
Kalman  filter  presented  in  Chapter  II.  The  user  simply  spec¬ 
ifies  the  dynamics  and  measurement  relations  for  his  system. 
SOFE  propagates  the  state  and  covariance  estimates  forward 
from  the  specified  initial  time,  using  a  fifth  order  Kutta- 
Merson  integration  algorithm  (6).  Updates  of  state  and  co- 
variance  estimates  based  on  available  measurements  are  pro¬ 
vided  at  user-specified  intervals.  The  user  can  specify 
any  number  of  measurements  to  be  incorporated  at  a  given  up¬ 
date  time.  Use  of  SOFE  greatly  reduces  the  amount  of  com¬ 
puter  programming  necessary  in  developing  a  Kalman  filter  and 
allows  the  user  to  concentrate  on  the  finer  details  of  his 
particular  problem. 

The  fifth-order  Kutta-Merson  integrator  implemented  in 
SOFE  requires  a  step-size  no  greater  than  approximately  two 
milliseconds  for  the  chosen  integration  tolerances.  The  in¬ 
tegrator  uses  a  variable  step  size  to  automatically  maintain 
the  integration  error  below  a  specified  value.  The  user  can 
specify  a  fixed  step  size  mode  if  exterior  factors,  such  as 
a  high  measurement  rate,  cause  the  step  size  to  remain  small 
regardless  of  dynamics.  If,  in  order  to  handle  severe  dynam¬ 
ics,  the  integrator  reduces  its  step  size  to  a  minimum  spec¬ 
ified  value  without  satisfying  error  tolerances,  an  integra¬ 
tion  failure  occurs  and  the  program  stops.  The  default 
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parameters  for  the  integrator  are  variable  step  size,  error 
tolerance  of  .0001,  maximum  step  size  of  1.0  E  +  9,  minimum 
step  size  of  .0001,  and  initial  step  size  of  .01.  We  are 
unsure  of  exact  vehicle  dynamics  but  expect  rapid  changes  in 
acceleration  over  a  very  short  time  interval.  We  desire  to 
allow  the  integrator  to  automatically  adjust  its  step  size 
in  order  to  reduce  integration  error  and  avoid  "stepping-over” 
any  fluctuations  in  the  solution.  A  variable  step  size  also 
reduces  computer  time.  For  these  reasons,  we  interpolate 
the  accelerometer  data  to  .002  second  intervals  for  integra¬ 
tion  purposes. 

S0FE  implements  user-supplied  data  records  from  one  ex¬ 
ternal  tape  and  expects  the  same  number  of  records  every  time 
the  tape  is  read.  The  external  data  provided  for  the  rocket 
car  includes  accelerometer,  range,  azimuth,  and  elevation 
data.  This  data  is  interpolated  to  .002  second  intervals 
using  the  cubic  spline  interpolator  implemented  in  SOFE.  Thus, 
accelerometer  data  for  integration  purposes  is  available  every 
.002  seconds  while  the  specified  measurement  update  interval 
of  .05  seconds  insures  that  actual,  not  interpolated  values 
for  the  radar  measurements  are  used. 

Accelerometer  Calibration 

The  first  run  through  the  data  using  SOFE  is  made  to 
check  the  calibration  of  the  accelerometer  data.  The  seven 
state  extended  Kalman  filter  as  implemented  in  SOFE  is  run 
without  incorporating  any  radar  measurements.  Integration 
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of  accelerometer  data  provides  state  estimates  of  x  position 
and  velocity.  We  desire  to  compare  the  estimate  of  x  veloc¬ 
ity  from  SOFE  at  FIM  trap  entry  time  to  the  actual  value 
recorded  by  the  FIM  system.  The  FIM  officials  had  used  a 
small  radio  transmitter  to  record  the  time  of  vehicle  entry 
into  the  first  trap.  This  constant  frequency  transmitter  was 
designed  to  stop  transmission  while  the  car  was  in  the  first 
trap,  and  then  begin  transmission  upon  trap  exit.  Based  on 
frequency  data  from  this  transmitter,  AFFTC  engineers  cal¬ 
culated  trap  entry  time  as  18.65  seconds  from  the  start  of 
the  run  (3).  The  estimate  of  state  x^,  velocity  along  the 
x-axis,  based  on  integrating  the  accelerometer  data  is  checked 
against  the  trap  reading  at  this  time.  The  vehicle  trajectory 
generated  by  integrating  the  accelerometer  data  is  presented 
in  Figs.  4.1  through  4.8.  The  speed  as  measured  by  the  FIM 
speed  trap  at  18.65  seconds  is  666.234  mph  or  977.1432  ft/sec. 
The  state  estimate  of  x-velocity  at  this  time  based  on  inte¬ 
gration  of  accelerometer  data  is  978.582  ft/sec.  The  FIM 
system  showed  that  the  vehicle  was  in  the  first  speed  trap 
for  -.108  seconds  (3).  At  an  average  speed  of  978.582  ft/ 
sec  over  the  trap  distance  of  105.6  ft,  we  confirm  a  time  of 
.1079  seconds  in  the  trap.  The  peak  speed  based  on  integra¬ 
tion  of  accelerometer  data  alone  is  1080.05  ft/sec  or  736.4 
mph.  This  velocity  occurs  at  16.85  seconds  into  the  run  and 
results  in  a  Mach  number  of  1.006  when  referenced  to  the  speed 
of  sound  of  1073.536  ft/sec.  The  accelerometer  data  indicates 
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the  rocket  car  was  above  the  reference  Mach  number  of  one 
for  approximately  1.25  seconds. 

Figure  4.2(a)  indicates  that  there  is  no  deviation  from 
center-line  due  to  our  model  for  y  velocity.  Without  any 
measurements,  no  estimate  of  y  position  or  velocity  can  be 
made.  Note  also  that  the  standard  deviations  of  state  esti¬ 
mates  increase  over  the  time  interval.  No  reduction  in  our 
initial  values  for  standard  deviation  is  possible  without  in¬ 
corporating  measurements . 

Calculation  of  Measurement  Noise  Variance 

Before  incorporating  the  available  radar  measurements  of 
range  and  azimuth  it  is  necessary  to  determine  the  errors  in¬ 
herent  in  these  observations.  Although  the  specifications 
for  the  radar  are  available,  one  has  sufficient  reason  to 
doubt  the  validity  of  these  numbers.  The  main  reason  for  con¬ 
cern  is  the  way  the  radar  is  used  to  track  the  vehicle.  This 
particular  radar  is  located  up  to  4.5  miles  from  the  vehicle 
on  a  hill  overlooking  the  lake  bed.  A  tracking  radar  normally 
used  to  track  airborne  vehicles  equipped  with  transponders 
is  being  used  in  a  "look-down"  mode  into  ground  clutter  at  a 
target  not  equipped  with  a  transponder.  As  described  in 
Chapter  I,  we  know  the  range  measurement  has  a  two-second-long 
interval  when  the  radar  picked  up  a  larger  vehicle.  In  addi¬ 
tion  to  the  known  error  in  the  range  measurement  this  data 
is  highly  suspect  for  the  reasons  outlined  above.  Certainly 
we  cannot  rely  on  the  specified  range  RMS  error  of  15  feet. 
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Nor  is  the  azimuth  measurement  expected  to  maintain  the 
specified  RMS  error  of  .005  radian. 

The  manual  tracking  adjustment  used  to  align  the  radar 
dish  with  the  car  is  subject  to  operator  errors.  An  AFFTC 
review  of  the  video-tape  from  the  television  monitor  used  to 
adjust  azimuth  tracking  rate  shows  that  the  operator  was  able 
to  keep  the  cross-hairs  on  the  vehicle  for  most  of  the  run 
(3).  However,  at  the  beginning  of  the  run  when  the  car  is 
accelerating  the  most,  and  immediately  after  engine  burn-out 
when  the  car  is  reaching  maximum  deceleration,  the  operator 
is  off  the  vehicle  by  up  to  three  car  lengths  (3).  This  very 
subjective  estimate  of  azimuth  deviation  corresponds  to  ap¬ 
proximately  108  feet.  At  the  minimum  range  of  18000  feet, 
this  deviation  contributes  an  error  in  azimuth  of  up  to  about 
.006  radian.  This  error  is  therefore  the  best  we  could  hope 
for,  assuming  no  inherent  azimuth  bias  error.  To  estimate 
the  accuracy  of  range  and  azimuth  measurements,  some  compari¬ 
son  of  actual  to  expected  values  is  necessary. 

Using  equations  (3-22)  and  (3-23)  for  radar  range  and 
azimuth  based  on  estimates  of  x  and  y  position,  we  desire  to 
compare  actual  measurements  of  range  and  azimuth  to  filter 
estimated  values  without  incorporating  radar  measurements. 

The  dynamics  model  using  accelerometer  data  is  used  to 
derive  the  estimated  range  and  azimuth  trajectories.  Based 
on  this  model,  y  position  remains  at  zero  while  the  estimate 
of  x  position  is  provided  by  twice  integrating  accelerometer 
output.  By  comparing  actual  to  estimated  measurements  we  not 


54 


only  can  find  the  noise  strength  of  the  measurements,  but 
check  our  equations  for  range  and  azimuth  used  in  the  filter. 

Figures  4.9(a)  and  4.9(b)  are  plots  of  actual  range  and 
estimated  range.  Figure  4.9(a)  includes  a  plot  of  the  dif¬ 
ference  between  these  two  values  for  radar  range.  It  is 
obvious  from  these  plots  that  between  16  and  18  seconds  the 
range  is  tracking  a  larger  vehicle  beyond  the  rocket  car.  It 
is  also  apparent  from  these  plots  that  the  range  measurement 
is  indeed  extremely  "noisy"  and  has  significant  errors. 

The  azimuth  measurement,  however,  appears  to  be  much 
better.  Referring  to  Fig.  4.10,  we  can  see  that  the  actual 
and  estimated  azimuth  values  are  very  close.  This  confirms 
our  assumption  that  the  radar  operator  did  a  good  job  of  track 
ing  the  car  in  azimuth. 

The  purpose  of  comparing  actual  to  estimated  measurements 
is  to  determine  realistic  values  for  the  diagonal  terms  in 
the  measurement  noise  matrix,  R(t^)  .  To  accomplish  this,  we 
sum  the  "residuals",  or  difference  between  actual  and  esti¬ 
mated  measurements,  over  the  entire  time  interval.  We  then 
calculate  a  mean  and  variance  for  the  residual  values  using 
the  following  equations 

N 

Mean,  u  =  1/N-l  S  r.  (4-1) 

res  i=1  i 

y  N  ?  ? 

Variance,  cT  =  1/N-l  E  (r.  -p  )  (4-2) 

res  i_1  i 

where  r^  is  the  residual  measurement  at  a  given  sample  time 
and  N  is  the  number  of  sample  periods  used  in  the  calculations 
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An  inherent  assumption  in  this  residual  analysis  is  that  we 

can  employ  the  principle  of  "ergodicity"  (4) .  This  principle 

says  simply  that  we  can  determine  the  statistics  of  a  random 

process  using  only  the  results  of  one  run  through  the  data 

rather  than  a  Monte  Carlo  analysis.  We  assume  that  the  sta- 
* 

tistics  of  the  residual  measurements  can  be  adequately  de¬ 
scribed  using  data  over  the  24  second  time  interval  of  in¬ 
terest.  The  residual  sequence  can  be  shown  to  be  a  white 
Gaussian  sequence  with  zero  mean  and  covariance 

c2res  =  [H(ti)P(ti’)HT(ti)+R(ti)]  (4-3) 


Thus  we  can  solve  for  an  estimate  of  the  measurement  noise, 

A 


R(t^),  for  each  observation  using: 


A 

R(ti) 


(4-4) 


Once  the  filter  calculated  values  for  error  variance  P(t^~), 

reach  "steady  state"  conditions,  the  second  term  in  (4-4) 

becomes  negligible  when  compared  to  the  residual  variance 
2 

o  res .  Thus,  we  use  the  calculated  residual  variance  from 
(4-2)  to  yield  the  initial  estimate  of  measurement  noise, 

A 

*{tL) . 

In  calculating  the  range  error  variance  no  residuals  in 
excess  of  300  feet  are  used.  This  effectively  "blocks"  the 
erroneous  range  measurements  between  16-18  seconds,  resulting 
in  the  use  of  450  of  487  total  measurements  between  0  to 
24.3  seconds.  All  of  the  azimuth  measurements  between  0  to 


24.3  seconds  are  used  to  calculate  the  mean  and  variance  of 
the  azimuth  residual. 

The  result  of  this  residual  analysis  indicates  an  error 

2 

variance  on  the  range  residual  of  22455.82  ft  .  The  azimuth 

.  -4 

residual  analysis  shows  an  error  variance  of  .5389  x  10 
2 

rad  .  These  variances  yield  calculated  RMS  errors  of  149.85 

_2 

ft  and  0.734  x  10  radian  for  the  range  and  azimuth  measure¬ 
ments,  respectively.  As  expected,  the  range  measurement 
error  is  much  greater  than  the  specified  error  due  to  the 
way  the  radar  is  used  to  track  the  vehicle.  The  azimuth 
error  is  very  close  to  the  first  guess  of  .006  radian. 

The  actual  and  estimated  radar  range  measurements  appear 
to  diverge  after  approximately  20  seconds.  Referring  to  Fig. 
4.9(b),  the  actual  radar  range  looks  fairly  good  between  20 
and  24  seconds.  This  divergence  of  filter  computed  and  ac¬ 
tual  range  values  caused  some  concern.  In  fact,  the  initial 
extended  Kalman  filter  runs  indicated  that  the  state  X£,  y 
position,  grew  unrealistically  to  approximately  300  feet  by 
24  seconds.  This  growth  in  X2  is  caused  by  the  range  measure¬ 
ment  which  appears  accurate  between  20  to  24  seconds.  It 
would  seem  that  perhaps  our  initial  assumption  of  the  test 
track  heading  is  incorrect.  In  fact,  if  the  actual  track 
heading  is  179  degrees  true,  not  due  south,  this  one  degree 
deviation  would  cause  an  approximate  change  in  y  position  of 
300  feet  if  the  car  is  17000  feet  down-track.  Thus,  it  is 
necessary  to  correct  the  relations  between  radar  measured 
components,  xf  and  y r ,  to  the  earth-fixed  coordinate  system. 
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Referring  to  Fig.  4.10(a),  the  coordinate  translation  be¬ 


comes  : 

xr  =  DELX  +  (cos  179° )x^  +  (sin  179° )x2 

yr  =  DELY  +  (sin  l°)x^  +  (cos  l°)x2 
or 

xr  =  DELX  -  .99985x1  +  .01745x2 
yr  =  DELY  +  .01745x1  +  .99985x2 


(4-5) 

(4-6) 


When  these  corrections  are  applied  to  (3-22)  and  (3-23)  and 
the  residual  plotting  and  variance  calculations  are  made, 
the  actual  and  estimated  range  measurements  are  much  closer. 
The  corrected  plots  of  range  and  azimuth  are  shown  in  Fig. 
4.11(a),  (b)  and  4.12.  With  these  corrections,  calculated 

residual  variance  for  range  becomes  16807.66  ft  and  azimuth 

•  ■  -4  2 

variance  reduces  to  .3573  x  10  rad  .  These  lower  variances 

result  in  RMS  errors  for  range  and  azimuth  of  129.88  feet  and 

.005977  radian,  respectively.  Note  that  we  have  calculated 

an  error  for  the  azimuth  which  matches  our  initial  guess  of 

.006  radian.  It  is  apparent  that  the  one  degree  correction 

is  closer  to  the  true  track  heading  and  reduces  the  majority 

of  modeling  error.  Thus,  we  have  computed  the  estimated 

A 

measurement  noise  strength  matrix,  R(t^): 


A 

R(ti) 


16807.66  0 

0  . 3573E-4 


(4-7) 
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ACTUAL, ESTIMATED,  AND  RESIDUAL  RANGE  MEASUREMENTS 
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Fig.  4.11(a) ,(b)  Residual  range  analysis  based  on  assumed 
track  heading  of  179  degrees  true. 


RflOfW  AZIMUTH  (RRDIRNS) 


ACTUAL  AND  ESTIMATED  RADAR  AZIMUTH 


Fig.  4.12  Residual  azimuth  analysis  based  on  assumed 
track  heading  of  179  degrees  true. 


Check  of  Filter  Implementation 

Based  on  the  vehicle  trajectory  generated  from  the  ac¬ 
celerometer  data,  we  desire  to  check  the  calculation  of  the 
observation  matrix,  H[t^ «x( t^_U .  The  model  we  are  using  to 
describe  the  propagation  of  x  velocity,  x^,  relies  heavily 
on  the  accelerometer  measurements  of  specific  force.  The 
filter  receives  very  good  information  on  the  behavior  of  x 
velocity  and  position,  but  relatively  poor  information  on  y 
velocity  and  position.  To  check  our  implementation  of  range 
and  azimuth  measurements  in  the  filter,  we  desire  to  remove 
accelerometer  specific  force  from  the  model  for  x^.  If  the 
range  and  azimuth  measurements  independently  track  the  basic 
trajectory  of  the  vehicle,  we  can  be  reasonably  sure  the 
observation  matrix  linearization  has  been  calculated  cor¬ 
rectly.  Note  that  we  are  looking  for  trends  and  not  specific 
confirmation  of  state  estimates.  Thus,  we  expect  the  range 
and  azimuth  to  show  a  similar  trajectory  for  the  car  to  that 
generated  by  the  accelerometer,  but  not  a  one-for-one  com¬ 
parison.  The  purpose  of  such  an  analysis  is  simply  to  in¬ 
sure  that  we  have  correctly  incorporated  the  measurements 
into  the  extended  Kalman  filter  state  and  covariance  update 
relations.  The  model  for  the  radar  range  and  azimuth  measure¬ 
ments  is  repeated  here  for  reference « 

zx  =  [(DELX-.99985x1+.01745x2)2+(DELY+.01745x1+.99985x2)2 

+zr2]i5+x6+v1 
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z2  =  arctan[(DELY+.01745x1+.99985x2)/(DELX-.99985x1 

+ .01745x2) ]+x7+v2 

Applying  (2-15)  to  these  equations  we  form  the  observation 
sensitivity  matrix,  HCt^jxft^-)]  .  it  is  the  model  for 
measurement  incorporation,  z(t^),  and  the  calculation  of 
HCt^;x(t^-)]  that  we  desire  to  check.  By  comparing  state 
trajectories  of  position  and  velocity  independently  obtained 
from  each  measurement,  we  hope  to  confirm  our  calculations 
and  modeling  techniques. 

The  primary  reason  for  deviations  between  measurement 
trajectories  and  the  accelerometer  profile  is  due  to  the 
rather  crude  model  we  substitute  for  the  x-velocity  state, 
x^ . .  We  now  choose  this  state  to  be  modeled  as  a  random  walk 
of  the  form 

x3  =  w3(t) 

where  the  driving  noise,  w3(t),  has  a  relatively  high  strength 
to  account  for  our  uncertainty  in  such  a  model  and  allow 
closer  tracking  of  actual  data.  Modeling  this  state  as  we 
have,  the  x  velocity  is  considered  a  constant.  Any  change 
to  its  value  based  on  measurement  updates  is  done  in  a  step¬ 
like  manner.  Therefore,  we  do  not  expect  to  get  exact  agree¬ 
ment  with  accelerometer  results . 

The  model  for  x-velocity  is  incorporated  into  a  six 
state  extended  Kalman  filter  where  accelerometer  error,  x^, 
has  been  removed.  The  filter  is  run  using  range  measurements 
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only,  azimuth  measurements  only,  and  both  measurements  com¬ 
bined.  The  results  of  these  three  runs  of  the  extended  Kalman 
filter  are  shown  in  Fig.  4.13  to  4.16.  Plots  are  made  of  the 
estimates  of  x  and  y  position  and  velocity.  Part  (a)  of  each 
figure  is  the  depicted  state  estimate  generated  using  range 
measurements  only.  Part  (b)  is  the  state  estimate  generated 
using  azimuth  measurements  only,  while  part  (c)  of  each  fig¬ 
ure  shows  the  result  of  combining  both  range  and  azimuth.  We 
see,  in  fact,  that  position  and  velocity  along  the  x-axis  be¬ 
have  as  we  would  expect  from  the  accelerometer  trajectory 
shown  in  Figs.  4.1  through  4.4.  However,  the  geometry  of  the 
radar  position  to  the  vehicle  is  not  conducive  to  accurate 
estimates  of  deviations  along  the  earth-fixed  y-axis.  The 
range  measurement  is  the  only  means  by  which  we  can  hope  to 
estimate  position  and  velocity  along  the  y-axis.  The  state 
estimates  of  these  values  are  subject  to  any  errors  in  the 
range  measurement  and  indicate  only  weak  observability  of  these 
states.  The  range  measurement  is  ignored  between  16-18  seconds 
to  account  for  the  known  error  in  this  measurement  during  this 
interval .  Due  to  the  geometry  of  the  problem,  the  range 
measurement  is  even  less  likely  to  track  the  vehicle  along 
the  x-axis  correctly.  This  can  be  seen  from  comparing  the 
plots  of  x-position  and  velocity  for  range  only  to  the  plots 
of  these  states  with  azimuth  only  and  both  measurements  com¬ 
bined.  From  these  plots,  it  is  apparent  that  the  azimuth 
does  a  credible  job  of  tracking  changes  in  position  and 
velocity  along  the  x-axis.  Conversely,  tl  azimuth 

66 


1 


o 

*  a 


DISTANCE  FROM  START  ALONG  X-AXIS 


fig.  4.13(a)  Random  walk  model  for  x^  -  Range  only 
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DISTANCE  FROM  START  ALONG  X-AXIS 


Fig.  4.13(b)  Random  walk  model  for  x^  -  Azimuth  only 
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DISTANCE  FROM  START  ALONG  X-AXIS 


4.13(c)  Random  walk  model  for  -  Range  and  Azimuth 
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Fig.  4.14(a)  Random  walk  model  for  -  Range  only 
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Fig.  4.14(c)  Random  walk  model  for  -  Range  and  Azimuth 
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VELOCITY  COMPONENT  IN  X-DIRECTION 


Pig.  4.15(c)  Random  walk  model  for  -  Range  and  Azimuth 
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Fig.  4.16(a)  Random  walk  model  for  -  Range  only 
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Fig.  4.16(b)  Random  walk  model  for  -  Azimuth  only 
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Fig.  4.16(c)  Random  walk  for  -  Range  and  Azimuth 
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measurement  can  tell  very  little  about  y-axis  changes  of 
position  and  velocity.  Again,  the  geometry  of  the  radar  and 
vehicle  track  causes  this  observability  problem. 

Figures  4.13  (c)  through  4.16  (c)  are  generated  by 
utilizing  both  azimuth  and  range  measurements.  Again,  the 
estimates  of  x  position  and  velocity  are  fairly  consistent 
with  the  estimates  of  these  values  generated  by  using  ac¬ 
celerometer  data.  Deviations  from  track  centerline  are  less 
severe  but  still  subject  to  range  measurement  errors.  It 
is  apparent  from  this  analysis  that  the  basic  trajectory  of 
the  vehicle  generated  using  accelerometer  data  is  reconfirmed. 
Thus,  we  conclude  that  the  extended  Kalman  filter  is  cor- 

r  A  -  1 

rectly  calculating  the  observation  matrix,  H(_t^;x(t^  )J, 
and  that  our  model  for  z(t^)  is  correct. 

Extended  Kalman  Filter  Performance  Analysis 

Having  calculated  measurement  noise  strengths  for  range 
and  azimuth  and  also  having  checked  the  calculation  of  the 
observation  matrix,  we  are  now  ready  to  run  the  extended 
Kalman  filter  modeled  in  Chapter  III.  The  initial  condition 
for  the  propagation  of  the  state  vector  is  chosen  as  zero- 
mean  Gaussian  random  variable  with  mean 

E  [x(tQ)>  ^  =  0  (4-8) 

and  covariance 

EC(x(t0)-x0)  (x(t0)-x0)T]  =  (4-9) 

where  the  initial  covariance  matrix  from  (3-16)  is  specified 
as  t 
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10000 

10000 

100 

100 

.0259 

225 


. 25E-6 


(4-10) 


The  initial  values  for  the  diagonal  entries  of  the  system 
noise  matrix,  Q.(t),  are  repeated  here  for  reference: 


au) 


0 

0 

0 

.01 
.0518 
1.0 
.  IE-7 


(4-11) 


The  extended  Kalman  filter  as  implemented  in  SOFE  is 
run  with  several  different  combinations  of  2(fc)>  and 

R(t^).  This  analysis  is  necessary  to  monitor  the  behavior 
of  the  error  variance  of  each  state  estimate  in  order  to 
check  filter  performance.  We  desire  to  obtain  the  lowest 
possible  error  variance  on  each  state  after  24  seconds  of 
filter  operation,  using  realistic  values  of  initial  covari¬ 
ance,  system  noise,  2(fc)»  and  measurement  noise,  R(t^). 

This  performance  analysis  can  be  thought  of  as  "tuning"  the 
extended  Kalman  filter.  The  results  obtained  from  a  large 
number  of  noise  and  initial  covariance  combinations  are  sum¬ 
marized  here. 
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The  initial  value  for  the  error  variance  matrix,  P, 
proved  to  be  adequate  for  all  the  states  except  range  and 

azimuth  biases.  The  initial  covariance  for  range  bias  and 

•  2  2 
azimuth  bias  had  to  be  increased  to  1600  ft  and  .5E-4  rad  , 

respectively,  to  see  some  reduction  in  the  error  variances 

of  these  states.  After  24  seconds  of  filter  operation,  the 

filter  only  reduces  the  error  variance  on  these  states  to 
2  2 

1190  ft  and  .88E-5  rad  ,  respectively  due  to  limited  obser¬ 
vability.  The  filter's  ability  to  estimate  accelerometer 
error  is  relatively  unaffected  by  changes  to  driving  noise 
on  this  state.  The  error  variance  on  the  filter  estimate  of 
accelerometer  error  coverges  very  quickly  to  its  minimum 
value  regardless  of  changes  in  driving  noise  or  initial  co- 
variance.  Likewise,  changes  to  the  pseudonoise  values  on 
states  x&  and  x7,  range  and  azimuth  bias,  do  not  effect  the 
filter's  ability  to  estimate  these  states. 

The  greatest  reduction  in  error  variance  on  all  the 
states  of  interest  results  from  the  recalculation  of  the 
measurement  noise  matrix,  R(t^),  based  on  analysis  of  the 
range  and  azimuth  measurement  residuals.  This  is  expected 
since  the  filter  now  "puts  more  stock”  in  the  measurements 
of  range  and  azimuth  than  it  did  for  higher  values  of  R(t^). 
Since  we  are  most  concerned  with  estimates  of  vehicle  veloc¬ 
ity,  x^  and  x^,  we  present  the  results  of  the  tuning  process 
on  the  error  variance  of  these  state  estimates.  The  minimum 
error  variances  for  these  states  occur  after  24  seconds  of 
filter  operation.  The  error  variance  of  x  and  y  velocity 
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estimates  and  associated  driving  noise  strengths  and  measure¬ 
ment  noise  is  presented  in  Table  I.  The  initial  variance 

2 

on  each  state  is  100  (ft/sec)  ,  or  ±10  ft/sec. 

TABLE  I 

Filter  Performance  Analysis 


Driving  Noise  q  Measurement  Noise,  R(t^)  P(t  =  24  sec) 


q4 

= 

.  2E-2 

22456/. 5388E-4 

P3 

— 

2.4 

q5 

= 

.0518 

•I 

II 

P4 

= 

6.36 

q6 

= 

10.0 

If 

II 

q7 

= 

.  5E-5 

II 

•  I 

q4 

= 

.  IE-2 

22456/. 5388E-4 

p3 

= 

3.193 

q5 

= 

.0515 

•I 

•I 

P4 

= 

6.43 

q6 

= 

1.0 

H  II 

q7 

.  IE-7 

II 

II 

q4 

= 

.  IE-3 

22456/. 5388E-4 

P3 

= 

2.57 

q5 

= 

.  3E-1 

•1  II 

P4 

= 

6.24 

q6 

= 

.10 

•I  II 

q7 

= 

.  IE-7 

•I  II 

q4 

= 

.  IE-7 

16808/. 3573E-4 

P3 

= 

2.27 

q5 

= 

.  3E-1 

II 

•  I 

p4 

= 

5.92 

q6 

= 

.01 

•I  II 

q7 

= 

.  IE-7 

•I  II 

The  combination  of  system  noise,  initial  covariance,  and 
measurement  noise  chosen  as  a  result  of  this  performance 
analysis  is  presented  here: 


78 


1 


Results  of  One  Iteration  of  the  Extended  Kalman  Filter 

State  estimates  and  error  standard  deviation  in  these 
estimates  for  the  "tuned”  extended  Kalman  filter  are  pre¬ 
sented  in  Fig.  4.17  through  Fig.  4.25.  These  plots  are 
generated  using  the  values  for  P^,  Q.(t),  and  R(t^),  specified 
above.  Figure  4.24(a)  is  a  plot  of  the  range  measurement 
residual  bracketed  by  the  residual  standard  deviation.  The 
azimuth  measurement  residual  and  associated  residual  standard 
deviation  is  plotted  in  Fig.  4.24(b).  It  should  be  noted 
that  a  "residual  monitoring"  routine  has  been  included  in 
the  basic  software  for  SOFE.  This  routine  calculates  the 
residual  standard  deviation  at  each  sample  time  from: 

ares  =  CH(ti)P(ti~)HT(ti)+R(ti)]Js 

and  compares  the  residual  measurement  to  this  value  for 

cf  „ .  If  the  residual  measurement,  6z(t-),  is  greater  than 
res  —  1 

3a  .  the  measurement  is  ignored.  As  a  result  of  residual 
res 
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Fig.  4.19(a)t(b)  Extended  Kalman  filter  -  Range  and  Azimuth 


VELOCITY  COMPONENT  ALONG  Y-AXIS 


Tam — ii.a-,1"  ibb,  .m — 'hiwuictu 


STD  DEV  OP  Y- VELOCITY 


Fig.  4. 20(a), (b)  Extended  Kalman  filter  -  Range  and  Azimuth 
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Fig.  4. 21(a), (b)  Extended  Kalman  filter  -  Range  and  Azimuth 
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Pig.  4.22(a) ,(b)  Extended  Kalman  filter  -  Range  and  Azimuth 
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Fig.  4. 23(a), (b)  Extended.  Kalman  filter  -  Range  and  Azimuth 
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Fig.  4. 25(a), (b)  Extended  Kalman  filter  -  Range  and  Azimuth 


monitoring,  40  measurements  between  16-18  seconds,  and  three 
measurements  near  six  seconds,  are  ignored. 

Figure  4.25  (b)  is  a  "blown  up"  version  of  Fig.  4.25 
(a)  which  is  the  magnitude  of  the  velocity  vector  converted 
to  Mach  number . 

The  plots  of  standard  deviation  for  range  bias  and 
azimuth  bias  indicate  weak  covergence  of  the  standard  devia¬ 
tion  for  these  states.  As  mentioned  previously,  this  is  due 
to  limited  observability  of  these  states.  The  standard  de¬ 
viations  of  the  other  state  estimates  show  good  covergence 
and  indicate  that  even  the  "noisy"  range  measurement  provides 
some  information  on  state  values  which  can  be  used  to  improve 
state  estimation.  From  the  plot  of  range  residual  standard 
deviation  it  should  be  apparent  that  the  residual  monitoring 
routine  is  rejecting  measurements  of  range  which  have  a 
residual  value  greater  than  3a  ,  or  approximately  390  feet. 

IT  co 

The  straight  line  segment  between  16-18  seconds  and  the  "spike" 
at  six  seconds  in  Fig.  4.24  (a)  show  where  the  range  meas¬ 
urements  have  been  ignored . 

The  large  deviation  in  y  position  shown  in  Fig.  4.18  (a) 
at  approximately  10  seconds  is  due  to  range  measurement 
errors  which  are  significant  between  6-10  seconds.  The  fil¬ 
ter  weights  these  measurements  lightly  due  to  relatively 
high  measurement  noise,  but  does  not  totally  reject  them. 

Thus,  the  estimated  y  position  shows  unrealistic  values  for 
this  state  due  to  inaccuracies  in  the  range  measurement.  It 
appears  from  Fig.  4.11  (b)  that  during  this  four  second 
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interval  between  6-10  seconds,  the  range  is  indicating  short 
of  actual  vehicle  position. 

The  extended  Kalman  filter  indicates  a  peak  speed 
estimate  at  16.85  seconds  of 


3  max 


=  1089.49  ft/sec 


4max 


=  -11.3765  ft/sec 


Calculating  the  scalar  speed  of  the  vehicle  from 

V  =  (x32  +  i42)^ 


we  find  a  maximum  scalar  velocity  of 


V  „  =  1089.55  ft/sec 
max 


resulting  in  a  maximum  Mach  number  of 


MmaY  *  1.015 

iHcix 

when  referenced  to  the  speed  of  sound,  a, 

a  =  1073.536213  ft/sec 

The  extended  Kalman  filter  indicates  that  the  rocket  car  was 
above  the  reference  speed  of  sound  for  approximately  1.9 
seconds . 

A  run  of  the  extended  Kalman  filter  using  only  azimuth 
measurements  was  made  to  check  the  influence  of  incorporating 
the  inaccurate  range  measurement.  The  estimate  of  velocity 
from  this  run  is  much  closer  to  the  estimate  of  x  velocity, 
x^,  obtained  from  integrating  accelerometer  data  without  in¬ 
corporating  any  measurements.  This  reconfirms  the 


1 


observability  problems  previously  mentioned  associated  with 
estimating  Xg  from  range  measurements.  The  azimuth  only 
run  of  the  extended  Kalman  filter  indicates  the  following 
vehicle  performance: 

x3max  =  1080.26  (16.85  sec) 

A 

Xg  =  977.081  (18.65  sec,  trap  entry) 

Time  above  Mach  1  =  1.30  sec 

i 

The  azimuth  only  run,  however,  does  not  provide  very  good 
estimates  of  y  position  or  velocity  due  to  the  observability 
of  these  states  from  azimuth  measurements  alone.  It  should 
be  noted,  that  the  apparent  error  in  x  velocity  caused  by 
observability  problems  using  range  measurements  is  corrected 
by  using  the  fixed-interval  smoother  algorithm.  The  time 
history  of  state  estimates  and  error  variances  from  the  run 
of  the  extended  Kalman  filter  incorporating  both  measure¬ 
ments  is  stored  for  use  in  the  smoother  algorithm.  The  re¬ 
sults  of  the  smoother  analysis  will  be  presented  in  the  next 
chapter . 
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V.  Optimal  Smoother  Results 


Chapter  II  describes  the  Meditch  form  of  the  optimal 
"fixed-interval"  smoother  algorithm  in  detail.  A  simple 
FORTRAN  program  was  written  by  the  author  to  incorporate 
this  smoother  algorithm.  The  output  of  the  extended  Kalman 
filter,  detailed  in  the  previous  chapter,  is  stored  for  use 
by  the  smoother  program.  The  FORTRAN  code  used  in  the  smoother 
program  is  listed  for  reference  in  Appendix  B.  The  results 
obtained  from  this  program  are  presented  in  this  chapter. 

The  output  of  the  optimal  smoother  at  the  initial  time, 
tQ,  is  used  to  correct  the  initial  conditions  for  another 
iteration  of  the  forward  extended  Kalman  filter.  Such  a 
"forward-backward"  iteration  scheme  is  used  to  correct  model 
errors  and  initial  conditions  of  the  extended  Kalman  filter. 
After  each  iteration  of  the  forward-backward  estimator  a 
comparison  of  state  values  is  made.  When  the  difference  in 
state  estimates  from  one  iteration  to  the  next  is  less  than 
some  arbitrarily  specified  value,  €,  the  estimator  is  said 
to  have  "converged".  Since  our  main  area  of  concern  is  in 
the  estimate  of  vehicle  speed,  we  compare  the  peak  estimates 
of  velocity  for  each  run  of  the  smoother.  When  the  differ¬ 
ence  between  peak  speed  estimates  from  one  iteration  to  the 
next  is  less  than  2  ft/sec  and  the  standard  deviation  of 
this  estimate  allows  at  least  a  99  percent  confidence  that 
the  vehicle  exceeded  the  speed  of  sound,  we  stop  the  itera¬ 
tions.  The  latter  requirement  for  standard  deviation 
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becomes  the  driving  factor  due  to  time  limitations  and  suf¬ 
ficient  confidence  after  two  iterations. 

One  Iteration  of  Smoother 

The  plots  of  the  "smoothed"  state  estimates  after  one 
smoother  iteration  are  presented  in  Figs.  5.1  (a)  to  5.7  (a). 
The  variance  for  each  state  estimate  calculated  by  the 
smoother  are  presented  in  Figs.  5.1  (b)  to  5.7  (b) .  Note 
that  the  smoother  works  "as  advertised"  in  reducing  the  error 
variance  of  state  estimates  when  compared  to  the  standard 
deviation  plots  of  the  extended  Kalman  filter  in  Chapter  IV. 
Some  very  interesting  results  are  obtained  and  need  to  be 
discussed.  The  plot  of  y  position.  Fig.  5.2  (a)  indicates 
that  the  vehicle  track  is  indeed  a  straight  line.  However, 
the  figure  indicates  that  y  position  begins  approximately 
250  feet  east  of  the  assumed  origin  and  decreases  in  a 
linear-manner  to  60  feet  east  of  centerline.  This  indicates 
that  the  starting  position  of  the  vehicle  is  displaced  east 
of  the  earth-fixed  coordinate  frame  and  that  the  assumed 
track  heading  of  179  degrees  true  is  incorrect  by  approxi¬ 
mately  .209  degrees.  The  smoother  estimated  state  values 
and  variances  of  x  and  y  position  at  the  start  of  the  run, 
tQ,  are  shown  to  be 

^l(to/tf)  =  167-78  ft  Pil(to/tf)  =  3746-5  ffc2 

A  9 

x2(t0/tf)  =  255.85  ft  P22(to,/tf^  =  3649.5  ft 
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Fig.  5. 7(a), (b)  Smoother  estimate  and  variance  of  azimuth  bias 
error  after  one  iteration. 


24.3  seconds: 


and  at  the  final  time,  tf  = 

x1(tf/tf)  =  16771  ft  Pnltf/tf)  =  3186.3  ft2 

A  9 

x2(tf/tf)  =  61.21  ft  p22(tf/tf)  =  1428<0  ft 

Figure  5.2  (a)  indicates  that  the  vehicle  track  is  indeed  a 

straight  line  but  not  aligned  with  the  x-axis  we  have  chosen. 
The  estimate  of  y  position  at  the  final  time  is  used  to  re¬ 
correct  the  track  heading  of  179  degrees.  The  y  position  at 
t^.  is  61.21  feet.  This  deviation  in  position  indicates  that 
we  have  "over-corrected"  track  heading  previously  by: 

x9(t,./t,)  ,  61.21 

tan  - £ — — ]  =  tan-  [ - 3  =  *209  degree 

Xl(tf/tf)  16771 

The  corrected  track  heading  now  becomes  179.209  degrees.  We 
can  apply  this  corrected  heading  to  the  measurement  rela¬ 
tion  (4-5)  and  (4-6)  so  that  the  bracketed  terms  in  these 
equations  now  become: 

(DELX  -  .9999x1  +  ,0138x2) 

(DELY  +  .0138x1  +  .9999x2) 

By  adjusting  the  track  heading  by  this  amount  and  recalcu¬ 
lating  the  residual  variance,  as  in  Chapter  IV,  we  should  be 
able  to  reduce  the  amount  of  measurement  noise,  R(t^). 

These  adjustments  should  help  reconfirm  our  knowledge  of 
vehicle  trajectory  and  allow  the  next  forward-backward  itera¬ 
tions  to  converge  to  the  "true"  state  values. 
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In  addition  to  resetting  the  initial  covariance  matrix, 
PQ,  to  reflect  smoother  covariance  P(to/t^.),  we  also  correct 
the  initial  conditions  of  the  state  vector  and  adjust  the 
amount  of  driving  noise,  Q.(t).  We  should  note  here  that 
several  methods  are  available  to  adjust  Q.(t).  The  backward- 
recursive  smoother  can  be  used  to  generate  an  estimate  of 

A 

system  noise  at  each  sample  time,  Q.(t^/tf)  (4).  One  can  also 
"tune"  Q.(t)  in  an  "off-line"  manner  and  allow  the  system 
noise  to  vary  over  the  time  interval  of  interest.  If  one 
has  knowledge  of  the  time-varying  nature  of  a  particular 
state,  this  knowledge  can  be  used  to  adjust  the  strength  of 
driving  noise.  For  example,  we  might  desire  to  relate  the 
amount  of  drivirftj  noise  on  the  azimuth  bias  error  state  to 
incorporate  our  Knowledge  of  radar  operator  tracking  perform¬ 
ance  versus  vehicle  acceleration.  As  a  side  note,  we  should 
also  menuon  that  the  smoother  algorithm  can  be  used  to  gener 
ate  an  estimate  of  the  applied  Controls,  u(t./tf).  An  esti- 
mate  of  the  controls  applied  t q  the  system  at  any  time,  t^, 
is  not  the  concern  of  this  analysis.  We  also  choose  to 
iteratively  adjust  2_(t)  and  use  constant  noise  levels  over 
the  time  interval  of  interest.  Such  an  iterative  adjustment 

A 

to  2.(t)  provides  adequate  smoother  performance  and  simplifies 
the  algorithm. 

Figure  5.4  (a)  is  a  plot  of  the  y  component  of  veloc¬ 
ity  obtained  after  one  iteration  of  the  smoother  algorithm. 

It  is  apparent  from  this  plot  that  the  y  velocity  indeed  be¬ 
haves  as  a  constant  with  a  value  of  approximately  -8.0  ft/sec 
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This  constant  negative  velocity  is  caused  by  the  erroneous 
initial  conditions  we  used  in  the  extended  Kalman  filter. 

As  seen  in  Fig.  5.4  (a)  the  smoother  is  able  to  detect  small 
deviations  from  the  constant  y  velocity.  This  is  seen  as 
small  "bumps"  in  the  plot  of  y  velocity.  Thus,  our  model  of 
random  walk  with  small  driving  noise  is  reconfirmed  and 
will  be  used  for  the  next  forward -backward  iteration. 

The  smoother  estimated  value  for  accelerometer  error, 
XgCti/tf),  can  be  seen  in  Fig.  5.5  (a).  The  error  appears 
to  grow  with  time  and  seems  to  be  related  to  velocity  (i.e., 
the  higher  the  velocity,  the  more  error) .  The  error  does 
appear  to  be  time-correlated  and  does  not  behave  as  a  constant 
bias  error.  Thus,  the  time-correlated  model  for  accelero¬ 
meter  bias  appears  to  be  valid.  We  can,  however,  adjust  the 
driving  noise  on  the  propagation  of  this  state  by  using  the 
smoother  calculated  steady  state  variance.  Figure  5.5  (b) 
shows  a  constant  variance  of  0. 14937 (ft/sec  )  after  only  a 
very  short  transient  period.  We  can  adjust  the  driving  noise 
on  this  state  by  using  this  steady-state  variances 

q5  =  P5(«)2/T  =  .029874  ( ft/sec2 ) 2/sec 

The  initial  variance  value  for  this  state,  PQ^,  is  also  set 

2  2 

to  .014937  (ft/sec  )  to  insure  a  stationary  accelerometer 
error  state  process. 

Figure  5.6  (a)  is  a  plot  of  smoother  estimated  radar 
range  bias  error,  ^(t^/t^).  This  figure  indicates  that 
<r  this  state  indeed  behaves  as  a  constant  with  only  slight 
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-177.49  ft. 


variations  from  the  initial  estimate,  ^(tQ/t^.)  = 

Thus,  a  reduction  in  the  amount  of  driving  noise  on  this 
state  appears  valid.  We  now  reduce  the  strength  of  driving 
noise  on  this  state  to  indicate  more  confidence  in  its  be¬ 
havior  as  a  constant.  The  value  of  q^  for  the  next  itera¬ 
tion  of  the, extended  Kalman  filter-smoother  combination  is 
\ 

obtained  by  ’'tuning"  the  extended  Kalman  filter  in  a  sensi¬ 
tivity  analysis.  The  new  value  for  is  determined  to  be: 

q6  =  .0001  f t2/sec 

This  amount  of  driving  noise  is  two  orders  of  magnitude  less 
than  the  value  used  in  the  first  iteration.  The  new  initial 
variance  for  x^  becomes: 

Pq6  =  1186.1  ft2 

and  the  new  initial  condition  on  this  state  obtained  from 
the  smoother  becomes : 

xo6  =  -177.49  ft 

Referring  to  Fig.  5.7  (a)  we  see  that  the  azimuth  bias 

A 

error,  x^(t^/t^),  does  not  behave  entirely  as  a  constant. 

It  appears  that  during  the  first  five  seconds  of  the  run 
the  bias  error  is  greatest  and  reduces  to  a  minimum  value  as 
the  vehicle  achieves  peak  speed  (minimum  acceleration) .  As 
the  vehicle  begins  to  decelerate,  the  azimuth  bias  error 
again  grows  to  a  larger  value.  This  result  is  consistent 
with  our  knowledge  of  radar  operator  tracking  error.  The 
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amount  of  driving  noise  on  this  state  appears  adequate  to 
allow  the  filter  to  track  deviations  in  its  value.  Thus, 
q -j  is  left  unchanged  for  the  next  iteration  of  the  forward- 
backward  estimator.  We  adjust  the  initial  condition  and 
variance  of  this  state,  as  before: 

Po7  =  .91619  E-5  rad2 

A 

xq7  =  .0066051  rad 

The  extended  Kalman  filter  and  smoother  are  used  again 
with  new  initial  conditions  on  the  states  and  adjusted  ini¬ 
tial  variance.  The  amount  of  driving  noise  on  accelerometer 
error,  x^,  and  range  bias,  Xg,  are  also  adjusted  for  the  next 
iteration  of  the  forward-backward  estimator.  As  a  result  of 
the  first  iteration  of  the  smoother  we  make  the  following 
adjustments  to  our  model  for  measurements  and  initial  con¬ 
ditions  : 

1)  The  track  heading  is  corrected  to  179.209  degrees 

true. 

2)  Based  on  the  corrected  model  for  measurement  in¬ 
corporation,  we  recalculate  measurement  residual  variance  to 

A 

reduce  the  estimate  for  R(t^),  R(t^)  . 

3)  Initial  conditions  on  the  state  vector  are  corrected 

A 

to  reflect  the  smoother  calculation  of  x(tQ/tf)  such  that 

~ 167 .78 
255.88 
A  1.1896 

=  -8.0108 
-.003834 
-177.49 
-.0064 
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4)  The  initial  covariance  matrix,  P  ,  is  adjusted  to 


^o  = 


—o' 

reflect  the  smoother  calculated  value  of  P(tQ/tf): 

3746.5 
3649.5 

2.315  0 

5.91 

0  .015 

1186.1 

. 916E-5 

5)  Finally,  the  system  noise  matrix,  g(t)  is  corrected 
to  indicate  increased  confidence  in  our  model  for  range  bias, 
and  adjusted  steady-state  variance  on  the  accelerometer  error: 

0 


a(t) 


0 

o  o 

g  . ie-7 

.02987 
.  1E-3 
.  IE-7 


The  residual  variance  analysis  detailed  in  Chapter  IV 

A 

now  produces  an  estimated  measurement  noise  matrix,  R(t^), 
such  that: 


R(tj_)  = 


16235.4  0 

0  . 3900685E-4 


The  range  residual  variance  is  calculated  from  415  of  487 
total  measurements  which  have  a  residual  magnitude  less  than 
300  feet.  These  variance  values  result  in  calculated  RMS 
errors  for  the  range  and  azimuth  measurements  of  127 .42  ft 
and  .006245  radian,  respectively. 

Incorporating  updated  values  for  measurement  noise, 
initial  state  and  variance  conditions,  and  system  noise,  we 


106 


rerun  the  extended  Kalman  filter.  The  results  of  this  run 
are  shown  in  Appendix  A.  The  state  and  covariance  time 
histories  are  stored  for  use  in  a  second  iteration  of  the 
smoother  algorithm. 

Second  Iteration  of  Smoother 

The  second  iteration  of  the  optimal  smoother  algorithm 
provides  refined  state  estimates  as  shown  in  Figs.  5.8  (a) 
to  5.14  (a).  The  basic  behavior  and  values  of  these  states 
remain  unchanged  from  the  first  iteration  of  the  smoother. 
The  smoother  state  estimates  at  the  initial  time,  t  ,  are 
refined  from  those  obtained  in  the  first  iteration. 

To  illustrate  the  convergent  properties  of  the  smoother 
we  present  a  comparison  of  state  estimates  and  variances  be¬ 
tween  the  first  and  second  smoother  iterations  in  Table  II. 
This  table  includes  the  percentage  difference  between  the 
two  iterations  for  each  state  and  variance  value  and  the 
overall  percentage  change  between  the  second  and  first  itera 
tion.  We  choose  to  compare  these  values  at  the  time  of  peak 
vehicle  speed  at  16.85  seconds.  Table  II  indicates  good 
reduction  in  error  variance  for  all  the  states.  This  is  due 
to  the  improved  initial  conditions  supplied  to  the  extended 
Kalman  filter  after  the  first  run  of  the  smoother.  The 
smoother  is  able  to  reduce  error  variance  from  the  first  to 
second  iterations  due  to  improved  state  estimates  from  the 
forward  filter.  The  convergence  of  the  velocity  state  esti¬ 
mates  and  reduction  in  error  variance  on  these  states  is 
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ACCELEROMETER  ERROR  ( 


Fig.  5-15  Vehicle  trajectory  in  earth-fixed  coordinate  - 
frame  of  reference (x2  plotted  against  -  feet) 
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Comparison  of  First  and  Second  Smoother  Iterations 


sufficient  to  stop  the  iterations  at  two.  Another  iteration 

of  the  smoother  using  refined  initial  conditions  for  the  ex- 

a 

tended  Kalman  filter  based  on  x(tQ/tf)  and  P(t  /tf)  was  not 
made  due  to  time  limitations  and  sufficient  confidence  in 
our  estimate  of  vehicle  performance.  A  third  iteration  of 
the  smoother  would  provide  only  limited  refinement  of  state 
estimates  and  better  reduction  in  error  variance.  However, 
it  will  be  shown  in  Chapter  VI  that  only  two  iterations  meet 
our  confidence  level  requirements. 

Figures  5.8  (a)  and  5.9  (a)  show  the  starting  position 
of  the  vehicle  to  be  approximately  400  feet  east  and  200  feet 
north  of  the  assumed  coordinate  origin.  The  x  velocity  esti¬ 
mate  (Fig.  5.10  (a))  at  tQ  again  indicates  that  the  vehicle 
has  already  started  down  the  track  at  the  assumed  initial 
time.  The  smoother  estimate  of  x  velocity  at  tQ  is  2.72  ft/ 
sec.  Comparing  this  to  the  accelerometer  only  run,  the  point 
in  the  radar  data  chosen  as  tQ  appears  to  be  in  error  by  ap¬ 
proximately  one  sample  period  or  .05  second.  The  y  velocity 
estimate  again  shows  a  constant  value  of  approximately  -8.33 
ft/sec.  This  constant  negative  velocity  is  caused  by  the 
orientation  of  the  test  track  with  respect  to  the  chosen  co¬ 
ordinate  frame  of  reference.  To  illustrate  the  actual  vehicle 
track  in  the  assumed  frame  of  reference.  Fig.  5.15  shows 
smoother  estimated  y  position  plotted  against  x  position. 

This  plot  indicates  the  estimated  trajectory  of  the  vehicle 
in  the  earth-fixed  coordinate  frame  of  reference.  The 
trajectory  is  indeed  a  straight  line,  but  not  along  the 
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x-axis  we  have  used.  However,  an  exact  starting  position 
and  track  orientation  are  not  the  goals  of  this  analysis. 

Our  main  interest  is  to  obtain  a  good  estimate  of  peak 
vehicle  speed  along  the  track.  We  are  not  concerned  with 
where  this  peak  speed  occurs  but  more  with  the  value  and 
error  of  this  estimate. 

The  second  iteration  of  the  smoother  yields  a  maximum 
velocity  of  1082.028  ft/sec  at  16.85  seconds  from  chosen 
initial  time.  The  scalar  speed  estimate  at  18.65  seconds, 

FIM  trap  entry  time,  is  975.043  ft/sec.  The  reason  we  do 
nor  get  better  agreement  between  the  smoother  and  FIM  esti¬ 
mates  of  velocity  at  the  trap  is  due  to  the  time  skew  in  the 
radar  data  previously  discussed.  The  time  scale  we  have  used 
shows  trap  entry  between  18.60  to  18.70  seconds  but  we  are 
not  sure  exactly  where  trap  entry  occurs  in  this  interval . 

The  behavior  of  the  error  states  indicated  by  the  first 
iteration  of  the  smoother  is  reconfirmed  by  the  second  itera¬ 
tion.  The  radar  range  bias  error  behaves  very  much  as  a 
constant  with  only  slight  deviations  from  a  steady  value  of 
-255.68  ft.  The  smoother  estimate  of  range  error  is  plotted 
in  Fig.  5.13  (a).  Radar  azimuth  error  in  Fig.  5.14  (a)  is 
again  shown  to  be  "slowly-varying”  over  the  24  second  inter¬ 
val.  Our  initial  assumptions  about  the  radar  operator  azimuth 
tracking  error  are  again  reconfirmed.  The  error  in  azimuth 
starts  out  high  as  the  operator  lags  behind  the  vehicle  due 
to  rapid  acceleration,  decreases  as  the  operator  "catches  up" 
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to  the  car  near  peak  speed,  and  again  increases  upon  vehicle 
deceleration  as  the  operator  “jumps"  ahead  of  the  car. 

The  accelerometer  error  state  estimate  shown  in  Fig. 

5.12  (a)  indicates  that  this  state  varies  between  approximately 
2 

0  to  0.08  ft/sec  or  0  to  .0025  g's.  Maximum  error  in  the  ac¬ 
celerometer  occurs  at  approximately  ten  seconds  into  the  run 
and  the  error  decreases  to  near  zero  by  24  seconds .  Perhaps 
the  behavior  of  this  error  state  can  be  explained  by  referring 
to  Fig.  3.2  which  is  a  plot  of  raw  accelerometer  data  in  g's 
versus  time.  Figure  3.2  indicates  that  maximum  sustained 
g's  on  the  vehicle  occur  between  0  to  10  seconds  and  slowly 
decrease  from  that  time  on.  It  appears  from  Fig.  5.12  (a) 
that  th<-  accelerometer  error  is  a  function  of  the  time  of 
application  and  level  of  sustained  g's  on  the  vehicle.  This 
figure  indicates  a  time-correlated  behavior  of  the  accelero¬ 
meter.  Such  behavior  may  have  been  adequately  modeled  as  a 
random  walk.  One  way  to  model  this  behavior  might  be  to  re¬ 
late  the  amount  of  driving  noise  on  the  accelerometer  error 
state  to  the  level  of  acceleration  units  at  any  given  time. 
Thus,  q,_,  could  be  modeled  as  time-varying  for  use  in  a 
random  walk  model  of  accelerometer  error.  Certainly,  it  can 
be  argued  that  a  correlation  time  of  one  second  is  too  short 
from  the  behavior  of  the  accelerometer  error  shown  in  Fig. 

5.12  (a).  Nevertheless,  the  forward-backward  iterations 
have  provided  better  information  on  the  "true"  behavior  of 
the  states  of  interest.  Another  iteration  of  the  smoother 
could  be  made  with  updated  initial  conditions  and  perhaps  a 
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different  model  for  accelerometer  error,  but  two  iterations 
have  provided  sufficient  reduction  in  variance  values  for 
our  purposes . 

Figures  5.8  (b)  through  5.14  (b)  are  plots  of  smoother 
calculated  error  variances  for  each  state.  These  plots  show 
that  the  backward  filter  is  able  to  reduce  the  errors  in 
state  estimates  from  those  obtained  from  the  first  iteration 
of  the  smoother.  After  two  iterations  of  the  forward-backward 
estimator,  the  error  in  state  estimation  is  reduced  by  an 
average  of  45%  over  that  obtained  in  the  first  iteration  of 
the  smoother . 

The  results  of  this  second  iteration  of  the  smoother  are 
now  used  to  test  the  hypothesis  that  the  rocket  car  did,  in 
fact,  exceed  the  reference  speed  of  sound.  This  will  be  shown 
in  detail  in  the  next  chapter. 
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VI .  Hypothesis  Testing 


The  previous  chapter  presented  the  results  of  the  ex¬ 
tended  Kalman  filter  -  fixed  interval  smoother  estimation 
scheme.  The  resulting  state  estimates  and  error  covariance 
after  two  iterations  of  the  smoothing  method  will  now  be 
evaluated  to  yield  the  best  estimate  of  peak,  rocket  car 
speed  and  a  confidence  level  for  this  estimate.  Before  we 
can  analyze  a  hypothesis  test  of  the  peak  vehicle  speed,  it 
is  necessary  to  calculate  the  scalar  speed  estimate  standard 
deviation. 

Development  of  Scalar  Speed  Standard  Deviation 

The  values  for  x^  and  x^,  x  and  y  velocity,  are  given 

A  A 

in  terms  of  mean  values,  x^t^/t^.)  and  x4(t^/tf)  and  variances 
P33 ( ti/tf ^  and  p44^ti//tf^  and  the  covariance 

Under  our  assumptions  of  approximately  Gaussian  error  models, 
these  mean  and  variance  values  completely  describe  a  two- 
dimensional  Gaussian  probability  density  function  which  pro¬ 
pagates  forward  in  time  from  the  initial  to  final  time.  The 
state  estimates  of  x  and  y  velocity  provide  the  components 
of  a  two-dimensional  conditional  mean  vector,  m,  the  magnitude 
of  which  is  the  estimate  of  scalar  speed  at  any  time,  t^. 

This  mean  vector  in  the  x-y  plane,  shown  in  Fig.  6.1,  locates 
the  peak  of  the  density  function.  Surfaces  of  "constant  like¬ 
lihood"  (4)  are  generated  by  passing  planes  through  the  den¬ 
sity  function  parallel  to  the  x-y  plane.  These  surfaces  are 
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Fig.  6.1  Velocity  vector  and  one  sigma  ellipse  of 
constant  likelihood 


ellipses  parallel  to  the  x-y  plane  and  when  viewed  from  above, 
the  one  sigma  surface  of  constant  likelihood  appears  as  in 
Fig .  6.1. 

The  covariance  matrix  for  the  two-dimensional  velocity 
density  function  determines  the  spread  of  the  density  about 
m  as  in  the  figure.  This  covariance  matrix  also  determines 
the  angular  orientation  of  the  "principal  axes"  (4)  of  the 
ellipses  of  constant  likelihood. 


a33  P34C34G43 

P43a43a34  a442 


(6-1) 


where  p  is  the  correlation  factor.  If  there  were  no  correla¬ 
tion  between  state  estimates  of  x  and  y  velocity  (i.e., 

P34  =  P43  =  0  ) »  the  covariance  matrix  would  be  diagonal, 

and  the  principal  axes  of  the  ellipse  would  be  parallel  to 
the  x  and  y  axes.  The  magnitudes  of  the  semimajor  and  semi¬ 
minor  axes  of  the  one  sigma  ellipse  are  determined  from  the 
eigenvalues,  of  the  covariance  matrix 


33  ~  ^7 

(6-2) 

44  =  ^2 

(6-3) 

where  the  primed  notation  indicates  the  lengths  are  defined 
in  the  principal  axes  frame  of  reference.  We  desire  to  re¬ 
late  the  calculated  scalar  speed  estimate  from  the  optimal 
smoother : 

[v(ti/tf  )|  =  [(x3(ti/tf))2+(x4(ti/tf))2]'2  (6-4) 
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to  its  associated  standard  deviation.  To  do  this  we  must 


find  the  associated  representation  of  the  normalized  veloc¬ 
ity  vector  estimate  in  the  frame  of  the  principal  axes  of 
the  one  sigma  ellipse  shown  in  Fig.  6.1.  The  velocity  vec¬ 
tor  in  the  x-y  coordinate  frame  can  be  represented  by: 


v 


xy 


=  v  1  +v  1 
xx  y  y 


(6-5) 


The  same  vector  can  be  specified  using  the  unit  vectors, 
and  i-n  the  principal  axes  frame  of  reference: 


v 


ele2 


+ 


(6-6) 


From  the  smoother-generated  variances  and  covariances 
for  x  and  y  velocity  we  solve  for  the  associated  eigenvalues 
and  eigenvectors  of  this  matrix.  The  square  root  of  the 
eigenvalues  of  this  2-by-2  matrix  will  be  shown  to  be  the 
length  of  the  semiminor  and  semimajor  axis  of  the  one  sigma 
ellipse.  The  associated  eigenvectors  determine  the  angular 
orientation  of  the  principal  axes  of  the  one  sigma  ellipse. 

We  solve  for  eigenvalues  and  eigenvectors  in  the  fol¬ 
lowing  manner.  For  simplification,  we  define  the  covariance 
matrix  of  the  velocity  states  at  any  time,  t^,  as: 


IWV  = 


^33  P34~l 


lP43 


44 


J 


(6-7) 


Forming  the  matrix  [Xl_-p]  and  solving  for  the  determinant 
we  can  find  the  eigenvalues  of  £34* 


Xl-P  = 


X_P33  -P34 

"P43  X-P44 


=  X  -<P33+P44)X+<P33P44-P34P43)  (6‘8) 


This  "characteristic  polynominal"  is  set  equal  to  zero  and 
the  resulting  roots  are  the  eigenvalues,  X^  and  of  the 

velocity  covariance  matrix.  Noting  that  this  covariance 
matrix  is  symmetric  so  that! 


P34  P43 

we  define  the  characteristic  polynominal  as: 

X2  -(P33-P44>X-( P33P44-P342> 


(6-9) 


(6-10) 


The  resulting  eigenvalues  are 

X 1  ’  2  3  C  ( P33+P44 }  P33+P44 }  { P33P44-P342  >  ^/2  (6-U) 
Substituting  these  eigenvalues,  X^,  into  the  matrix  [ X .  I_-P ] 
and  noting  that  the  eigenvectors  are  in  the  "null  space"  of 
this  matrix,  so  that 


Xi_P33 


-P 


34 


-P„  X.-P, 


e . 

xl 

e  ■ 

.  2. 

=  0 


(6-12) 


results  in  two  equations  for  each  eigenvector  of  the  form: 

a  e.  +  be.  =  0  (6-13) 

n  12 

either  one  of  which  is  related  to  the  other  by  a  constant. 
Thus,  we  have  one  equation  and  two  unknowns  from  which  to 
solve  for  the  eigenvector,  ,  for  a  given  eigenvalue,  X^. 

We  need  another  relation  between  the  components  of  the  eigen¬ 
vector  in  order  to  solve  for  the  individual  elements.  Noting 
that  the  "normalized"  eigenvector  is  the  unit  vector  we  find 

the  other  equation:  2  ? 

ei  +  e.  z  =  1  (6-14) 

1  12 
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With  two  equations  and  two  unknowns  we  can  solve  for  the  in¬ 
dividual  components  of  each  eigenvector. 

The  resulting  eigenvectors  determine  the  angular  orienta¬ 
tion  of  the  principal  axes  of  the  one  sigma  ellipse.  The 
angle,  01#  between  the  x-axis  and  the  semiminor  axis  of  the 
ellipse  (associated  with  x  velocity  error  variance)  can  be 
found  from: 

=  tan"1  (e12/ell)  (6-15) 

where  e^  and  e^  are  the  components  of  the  eigenvector  which 
describes  the  orientation  of  the  semiminor  axis.  This  angle 
specifies  a  coordinate  transformation  matrix,  L,  for  a  rota¬ 
tion  about  the  z  axis,  such  that: 

cos  0  s in  0  0 

L  =  -sin©  cos©  0  (6-16) 

1  L  0  0  :- 

This  coordinate  transformation  matrix  is  used  to  relate  the 
orientation  of  the  normalized  velocity  vector  to  the  princi¬ 
pal  axes  of  the  one  sigma  ellipse.  Transforming  this  veloc¬ 
ity  unit  vector  into  the  frame  of  reference  of  the  principal 
axes,  we  obtain  a  length  from  ellipse  center  to  the  one  sigma 
ellipse  in  the  direction  of  the  velocity  vector.  It  is  the 
magnitude  of  this  length,  0vel  in  Fig.  6.1,  which  determines 
scalar  speed  standard  deviation. 

Calculation  of  Peak  Scalar  Speed  and  Standard  Deviation 

After  two  iterations  of  the  smoother  algorithm  we  arrive 
r  at  the  following  estimates  of  x  and  y  velocity  at  16.85 
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seconds.  As  in  the  previous  iteration,  this  time  is  found 
to  be  the  point  at  which  peak  x  velocity  is  obtained.  The 
associated  error  variances  and  covariances  for  x3  and  x4  at 
this  time  are  also  givens 


Velocity  State  Estimates 


Covariances 


x3(16.85/tf)  =  1081.996  fps  ?33  =  1.117148  (fps)' 
x4(16.85/tf)  =  -8.335287  fps  P44  =  2.140282  (fps)' 


P34  -  P43  .071  (fps) 


The  magnitude  of  the  velocity  vector  or  scalar  speed  estimate 


v  -  1082.028  fps 


(6-17) 


From  the  covariances  for  the  velocity  state  estimates  at 
16.85  seconds  we  form  the  two-dimensional  covariance  matrix: 


1.117148 


2.140282 


(6-18) 


The  eigenvalues  and  associated  eigenvectors  for  this  matrix 


are  found  from  (6-11),  (6-13),  and  (6-14): 


\l  =  1.11328 


el  = 


X2  =  2.14453125  e2  = 


99852 


.0544 


.06894 

.99762 


(6-19) 


The  lengths  of  the  semiminor  and  semimajor  axes  of  the  one 


sigma  ellipse  of  constant  likelihood  become: 


1 


(6-20) 


=  1.05512  ft/sec 
VI2  =  1.464422  ft/sec 

The  eigenvector  associated  with  X^  determines  the  angular 
orientation  of  the  semiminor  axis  of  the  one  sigma  ellipse 
in  the  x-y  coordinate  frame  of  reference: 

0  =  tan-1( .0544/. 99852)  =  3.11843  degrees  (6-21) 

This  angle,  0^,  defines  a  direction  cosine  matrix,  L,  in 
order  to  transform  the  nomalized  velocity  vector 


v  =  1081 . 996 i  -8 . 335287i  =^.  99997i  - .  0077i 
x  y  x  y 


(6-22) 


into  the  frame  of  reference  of  the  one  sigma  ellipse  such 


that 


.99852  .0544  0 

L  =  -.0544  .99852  0 


(6-23) 


i  =  .99852e1  +  .0544e2 


i  =  -.05446,  +  . 99852e0 
y  1  2 


(6-24) 


The  normalized  velocity  vector  in  the  principal  axes  frame 
becomes 

v-  -  =  . 99997 ( . 998525. + . 0544eo ) - . 0077 ( - . 0544e. + . 9985 2e7 ) 

ele2  12  12 

=  .9989i1+.0467e2  (6-25) 

The  angle  of  this  normalized  vector  with  respect  to  the  x-axis 
becomes : 


)2  =  tan-1 ( .0467/. 9989)  =  2.6767  degrees 


(6-26) 
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Subtracting  (6-26)  from  (6-21)  we  find  the  angle  between  the 
e^  eigenvector  and  the  velocity  vector  in  the  principal 


axes  f  rame  s 


-  9 2  ~  *44172  degrees 


(6-27) 


This  knowledge  will  become  helpful  in  a  moment. 

We  need  to  find  the  distance  from  the  one  sigma  ellipse 
center  to  the  ellipse  itself,  in  the  direction  of  the  veloc¬ 
ity  vector,  v-  -  .  To  accomplish  this,  we  have  the  general 
ele2 

equation  for  any  point  (x^,  x2)  on  the  ellipses 


2  2 

X1  ^2_ 

h  X2 


(6-28) 


where  c  is  a  constant.  For  the  one  sigma  ellipse  c  is  equal 
to  one.  We  also  have  the  familiar  relationship  between  two 


vectors  v^,  y2*. 


V1  *  v2 
rTZ — i - r=-|-  cos  Y 


(6-29) 


where  y  is  the  angle  between  the  vectors .  Let  x^  and  x2 
describe  the  coordinates  of  the  point  where  the  line  from 
ellipse  center  in  the  direction  of  the  transformed  velocity 
unit  vector  intersects  the  ellipse.  We  apply  (6-28)  and 
(6-29)  to  the  point  on  the  ellipse  described  by  these  points 
x^  and  x2*  In  (6-29)  we  are  interested  in  the  ij  eigenvector 
direction  of  length  and  the  transformed  velocity  vector 


V  =  Xjej  +  x2e2 


(6-30) 


i 


/ 2  2  . 

of  length  Vx^  +  x2  •  It  is  this  length  which  describes  the 
one  sigma  deviation  of  the  scalar  speed  estimate.  From 


(6-28) 


(6-31) 


and  from  (6-29) 


VX1e1  .  [x1e1+x2e2] 


cos  y  =  cos (.44 172) 


xlel«12e2 


27  I 


Vx1"+x2i-  =  .99997 


129 .095572X- 


From  (6-31) 


we  find 


and  from  (6-32) 


16665. 67x22  +  x22  =  1 


x2  =  .008173 


Thus , 


x1  =  1.055098 


x2  =  1.05513 


.99997 


(6-32) 


(6-33) 


(6-33) 


This  value  is  the  standard  deviation  of  the  scalar  speed  esti¬ 
mate,  1082.028  ft/sec.  Note  chat  we  expected  the  standard 
deviation  of  the  velocity  estimate  to  be  only  slightly  higher 
than  the  one  sigma  deviation  of  the  x  velocity  estimate  due 
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to  the  small  off-diagonal  terms  of  the  matrix,  and  rela¬ 
tive  magnitudes  of  the  x  and  y  velocities. 

Hypothesis  Test  of  Peak  Speed  and  Conf idence  Level 

We  are  now  prepared  to  apply  a  hypothesis  test  of  the 
peak  vehicle  speed.  We  wish  to  test  the  hypothesis  that  the 
peak  speed  estimate  is  above  the  reference  speed  of  sound,  a 

a  =  1073.536213  fps 

From  our  assumptions  of  Gaussian  models,  the  estimate 
of  peak  speed  and  associated  covariance  describe  a  condi¬ 
tional  normal  distribution.  The  appropriate  one-sided  "con¬ 
fidence  interval"  (8)  for  this  hypothesis  test  is  given  by 
(8): 

x  -  z  (confidence  level)  a(x)  >  u  (6-34) 

where 

x  -  mean  of  normal  distribution  =  1082.028 
z  -  area  under  the  standard  normal  distribution 
curve 

o(x)  -  standard  deviation  of  normal  distribution  = 
1.05513 

4  -  lower  bound  of  confidence  interval  =  1073.536213 
For  x  =  1082.028  fps,  4  =  1073.536213,  a(x)  =  1.05513: 

z  ( conf  idence)  =  =  8.048  (6-35) 

To  eight  significant  figures  (6-35)  yields  a  probability 
that  the  vehicle  was  below  Mach  one  of : 


1  -  P(x)  =  3.4346578  E-15  (6-36) 

from  which  P(x),  the  probability  that  the  vehicle  was  above 


Mach  one,  is  found  to  be: 

P(x)  =  .9999999999999965653422  (6-37) 

For  all  intents  and  purposes,  we  have  achieved  a  proba¬ 
bility  or  confidence  level  of  one  that  the  vehicle  exceeded 
the  reference  speed  of  sound.  Of  course,  this  confidence 
level  is  based  on  the  assumptions  and  modeling  techniques  we 
have  used  in  this  analysis.  Such  a  high  confidence  after 
only  two  iterations  of  the  forward-backward  estimator  illus¬ 
trates  the  power  of  optimal  smoothing  theory  in  post-run 
data  analysis. 


VII 


Conclusions  and  Recommendations 


Conclusions 

Two  iterations  of  the  forward-backward  estimator  have 
provided  improved  state  estimates  and  lower  error  variance 
than  is  possible  with  a  forward  filter  only.  The  backward 
recursive  fixed-interval  smoother  provides  updated  initial 
state  and  variance  values  which  yield  improved  state  esti¬ 
mates  from  the  forward  extended  Kalman  filter.  The  forward- 
backward  estimation  method  has  reduced  the  error  variance 
from  the  forward  filter  by  more  than  half  after  only  two 
iterations.  Based  on  the  second  iteration  of  the  smoothing 
algorithm  we  make  some  general  comments. 

The  improved  initial  time  conditions  from  the  smoother 
indicate  an  inaccurate  assumption  of  the  starting  position 
of  the  run.  We  used  an  origin  based  on  very  inaccurate  range 
measurements.  It  appears  that  the  initial  x  and  y  position 
are,  in  fact,  displaced  approximately  200  feet  north  and  400  ft 
east  of  the  assumed  starting  position.  The  vehicle  is  also 
already  moving  at  our  assumed  initial  time.  This  can  be  ex¬ 
plained  by  an  error  in  the  initial  time  chosen  for  the  radar 
measurements.  The  point  in  time  in  the  radar  data  chosen  as 
t  appears  to  be  off  by  approximately  0.05  second.  In  other 
words,  the  actual  starting  time  of  the  run  is  about  one  radar 
data  sample  before  the  time  chosen  as  the  starting  time. 

This  causes  the  smoother  to  estimate  an  off-zero  velocity  at 
our  declared  initial  time.  This  also  explains  why  we  do  not 
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see  exact  correlation  between  filter  estimated  velocity  and 
FIM  trap  speed  at  18.65  seconds.  On  the  time  scale  we  have 
used,  the  FIM  trap  occurs  somewhere  between  18.60  and  18.70 
seconds.  In  spite  of  this  time  skew,  the  estimated  values 
are  close  enough  to  trap  speed  to  allow  the  comparison  shown 
in  Table  III . 

The  actual  starting  position  of  the  vehicle  is  really 
not  the  information  we  desire.  We  set  out  primarily  to  get 
the  best  estimate  possible  of  peak  vehicle  speed  no  matter  at 
what  time  or  where  on  the  track  this  occurs .  In  terms  of 
vehicle  velocity,  the  smoothing  algorithm  used  in  this  analy¬ 
sis  after  two  iterations  has  provided  excellent  convergence 
to  the  "true"  peak  speed. 

We  also  have  come  close  to  the  maximum  velocity  estimate 
obtained  on  the  day  of  the  run  from  AFFTC  radar  data  analysis. 
It  appears  that  the  AFFTC  method  used  to  correct  erroneous 
range  data  was  valid  and  even  averaging  only  three  radar 
points  came  very  close  to  the  "true"  peak  speed.  We  now 
summarize  the  estimates  obtained  of  the  peak  vehicle  speed 
by  AFFTC,  accelerometer  data  only,  and  one  and  two  iterations 
of  the  forward-backward  smoothing  method  incorporating  range 
and  azimuth  measurements.  These  velocity  estimates  are  con¬ 
verted  to  Mach  number  using  the  calculated  reference  speed 
of  sound  of  1073.536213  ft/sec.  These  results  are  summarized 
in  Table  IV. 

Two  iterations  of  the  optimal  smoother  also  provide  some 
information  on  the  behavior  of  the  error  states  of  the 
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TABLE  III 


Speed  Estimate  at  FIM  Trap  (18.65  seconds) 


Method 

Feet/sec 

MPH 

FIM  Recorded  Speed 

977.1432 

666.234 

Accelerometer  Data  Only 

978.582 

667.215 

First  Run  of  Extended  Kalman 
Filter 

979.22 

667.649 

First  Run  of  Smoother 

972.91 

663.35 

Second  Run  of  Extended  Kalman 
Filter 

978.568 

667.206 

Second  Run  of  Smoother 

975.043 

664.802 

TABLE  IV 

Peak  Scalar  Speed 

Estimates 

Method  Used 

Feet/sec 

MPH 

MACH 

Time  Above 
Mach  1 

AFFTC  Computer  Analysis 
of  Corrected  Radar  Data 

1084.835 

739.66 

1.0105 

N/A 

Integration  of  Longitud¬ 
inal  Accelerometer  Data 

1080.05 

736.4 

1 .006 

1 . 25  sec 

First  Run  of  Extended 
Kalman  Filter 

1089.5 

742.84 

1.0149 

2.0  sec 

First  Run  of  Smoother 

1080.006 

736.34 

1.006 

1.25  sec 

Second  Run  of  Extended 
Kalman  Filter 

1086.71 

740 . 94 

1.0123 

1.8  sec 

Second  Run  of  Smoother 

1082.028 

737.75 

1.008 

1.4  sec 
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accelerometer  and  radar.  The  accelerometer  error  varies  be¬ 


■r 


tween  0  and  .08  ft/sec^(0  to  .0025  g's)  during  the  24  second 
time  interval  of  interest,  achieving  its  maximum  value  at  ap¬ 
proximately  ten  seconds  into  the  run.  From  ten  seconds  on, 
this  error  slowly  decreases  to  approximately  zero  by  24  sec¬ 
onds.  The  sustained  g's  on  the  vehicle  are  fairly  high 
(Fig.  3.2)  up  to  ten  seconds  and  then  begin  to  decrease  after 
this  time.  It  would  appear  the  accelerometer  error  is  a 
function  of  the  length  of  time  sustained  g's  are  applied  to 
the  accelerometer  and  the  magnitude  of  these  acceleration 
units.  Depending  on  one's  definition  of  "slowing-varying" , 
one  could  make  a  case  for  using  a  random  walk  model  for  the 
accelerometer  error.  Certainly,  it  could  be  argued  that  a 
correlation  time  of  one  second  is  too  short  for  the  behavior 
of  this  error .  Another  study  of  the  rocket  car  data  could 
use  on-line  "tuning"  of  the  system  noise  matrix,  Q.(t),  by 
allowing  the  smoother  to  calculate  an  estimate  of  its  value 

A 

over  time,  Q.(t^/tf).  One  could  also  calculate  smoother  esti- 

A 

mated  inputs  u(t^/t^.)  .  In  terms  of  the  rocket  car  analysis, 
such  estimation  of  accelerometer  input  at  any  time  t^  based 
on  the  entire  measurement  time  history  would  yield  improved 
state  estimation.  Nevertheless,  the  smoother  has  provided  a 
better  "glimpse"  of  the  "true"  behavior  of  this  state  than 
is  available  from  a  forward  filter  only,  especially  without 
more  knowledge  about  inherent  accelerometer  errors. 

The  random  walk  models  for  radar  range  and  azimuth  bias 
errors  prove  to  be  very  adequate.  These  errors  are  shown  as 
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nearly  constant  over  the  24  second  interval  with  only  slow 
changes  from  the  behavior  of  true  constants.  Certainly,  we 
are  somewhat  surprised  at  the  magnitude  of  the  radar  range 
bias  error.  Considering  the  size  of  the  vehicle  being  tracked, 
lack  of  transponder,  ground  clutter,  and  distance  from  the 
radar  site,  it  is  conceivable  that  the  radar  range  has  a 
large  inherent  error.  The  azimuth  bias  error  behaves  as  we 
expected  based  on  our  knowledge  of  operator  tracking  perform¬ 
ance.  Perhaps  "bias"  is  a  misnomer,  as  most  of  the  error  in 
azimuth  is  operator-induced.  The  azimuth  error  shows  that 
the  operator  lags  behind  the  vehicle  initially  but  is  able  to 
regain  good  tracking  as  the  acceleration  decreases.  After 
engine  "flame-out"  at  approximately  18  seconds,  the  azimuth 
error  again  increases,  indicating  the  operator  has  probably 
"jumped"  ahead  of  the  vehicle.  For  another  iteration  of  the 
estimator,  a  better  description  of  the  azimuth  error  could 
be  used.  One  could  relate  the  strength  of  driving  noise  on 
this  state,  q ~j,  to  the  acceleration  of  the  vehicle.  When  ac¬ 
celeration  is  high  q-,  would  be  increased.  The  amount  of  driv¬ 
ing  noise  would  decrease  as  vehicle  acceleration  decreases. 

This  analysis  has  shown  that  state  estimation  can  be 
significantly  improved  if  the  estimation  algorithm  has  access 
to  future  measurements.  This  is  the  real  benefit  of  a 
smoother  algorithm  in  post-run  data  analysis.  The  method 
used  in  this  analysis  requires  a  straightforward  incorpora¬ 
tion  of  existing  theory  and  available  software  with  only 
limited  additional  programming  required.  The  forward-backward 


iteration  scheme  is  a  simple  yet  effective  way  to  provide 
improved  state  estimates  in  an  "off-line"  application.  It 
can  be  applied  to  almost  any  system  of  interest  no  matter 
what  dimension,  with  only  computer  workload  becoming  a  driv¬ 
ing  factor.  Good  dynamical  and  error  models  are  a  require- 
ment,  but  not  a  necessity.  The  iterative  method  used  in 
this  analysis  can  help  "fine-tune"  very  simplified  models  to 
provide  improved  state  estimates. 

The  reader  familiar  with  estimation  of  unknown  para¬ 
meters  using  a  "maximum  likelihood"  estimation  technique  may 
wonder  if  such  a  technique  could  have  been  employed  in  this 
analysis.  The  answer  is  a  guarded  "yes"  if  we  can  make  some 
valid  assumptions.  The  inherent  assumption  in  maximum  likeli¬ 
hood  estimation  is  that  the  parameters  to  be  identified  can 
be  accurately  modeled  as  constants  over  some  time  interval 
of  interest.  In  this  analysis,  we  are  concerned  with  accur¬ 
ate  estimates  of  state  values  at  discrete  points  in  time  (i.e., 
peak  speed  at  some  time,  t^).  If  we  assume  the  parameters 
affecting  this  problem,  such  as  accelerometer  and  radar  errors, 
are  constant  over  time,  a  maximum  likelihood  estimation  algo¬ 
rithm  will  yield  a  best  fit  of  a  constant  to  the  data.  If 
the  parameters  are  not  true  constants,  a  better  (non-constant) 
model  would  inherently  allow  better  estimation  accuracy. 

This  analysis  has  shown  that  the  error  states  do  not  behave 
as  constants.  Therefore,  one  cannot  accurately  model  these 
errors  as  constant  unknown  parameters  for  implementation  in 
a  maximum  likelihood  estimation  algorithm  without  non-negligible 

138 


1 


estimation  performance  degradation.  Without  a  priori  informa¬ 
tion  on  the  behavior  of  the  errors  to  be  modeled,  the  deci¬ 
sion  was  made  to  model  the  errors  as  states  in  the  forward 
extended  Kalman  filter-backward  smoother  estimator. 

If  there  is  a  "bottom-line"  to  this  analysis  it  has  to 
do  with  the  peak  speed  of  the  Budweiser  Rocket  Car  on  17 
December  1979.  Rather  than  "eyeballing"  the  peak  speed  of 
the  car  based  on  poor  data,  we  are  able  to  provide  an  esti¬ 
mate  of  the  speed  and  a  confidence  level  for  our  estimate. 

In  fact,  after  only  two  iterations  of  the  forward-backward 
smoothing  technique,  we  can  state  with  probability  of  nearly 
one  that  the  vehicle  did  achieve  the  reference  speed  of 
sound,  based  on  the  assumptions  and  modeling  techniques  used 
in  this  analysis. 

Recommendations 

The  position  estimates  and  off-zero  velocities  calcul¬ 
ated  by  the  smoother  at  the  initial  time,  tQ,  indicate  a 
rather  poor  choice  of  origin  for  the  vehicle  frame  of  refer¬ 
ence.  Relying  on  the  radar  to  provide  a  good  initial  "fix” 
of  vehicle  position,  no  matter  how  long  the  radar  is  aimed  at 
the  vehicle,  is  only  "wishful  thinking".  Perhaps  a  better 
origin  could  have  been  located  at  remote  camera  site  A8 
shown  in  Fig.  1.1.  This  point  has  been  surveyed  and  "exact" 
latitude  and  longitude  coordinates  of  both  A8  and  the  radar 
site  are  known.  Using  these  coordinates  one  could  calculate 
a  much  better  DELX  and  DELY  from  which  to  reference  changing 
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radar  measurements  to  vehicle  motion  in  the  frame  of  refer¬ 
ence.  If  the  information  on  track,  alignment  with  respect  to 
A8  is  correct,  a  frame  of  reference  at  A8  should  indicate  a 
vehicle  track  parallel  to  the  x-axis.  The  starting  position 
of  the  vehicle  within  a  frame  of  reference  centered  at  A8  is 
still  unknown,  however,  and  only  iterative  methods  could 
"zero-in"  on  the  "true"  starting  position. 

An  adjustment  could  be  made  in  the  initial  time  chosen 
in  the  radar  data  to  find  the  "true"  sample  time  as  the 
vehicle  starts  to  move.  At  20  samples  per  second,  however, 
one  can  only  get  within  0.05  second  accuracy.  Also,  the 
radar  data  is  constant  until  the  azimuth  suddenly  increases 
very  rapidly.  We  chose  one  sample  time  before  the  first 
change  in  azimuth  as  the  initial  time.  One  could  "back-up" 
the  radar  data  until  the  smoother  estimate  of  x-velocity  at 
tQ  approaches  zero. 

Other  possibilities  for  further  study  include  some  off¬ 
line  tuning  of  system  noise  to  account  for  the  time-varying 
nature  of  accelerometer  and  azimuth  errors.  One  way  to  ac¬ 
complish  this  might  be  to  use  the  smoother  estimate  of  sys- 

A 

tern  noise,  Q/t^/t^),  based  on  the  measurement  data  to  provide 
a  time  history  of  driving  noise  for  each  of  the  affected 
states.  This  would  provide  the  forward  extended  Kalman  fil¬ 
ter  with  improved  knowledge  of  state  behavior . 

Finally,  it  might  be  beneficial  to  allow  the  smoother 
to  calculate  an  estimate  of  the  applied  controls,  in  this 
case  accelerometer  specific  force.  This  smoother  estimated 
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control  input,  utt^/t^.),  could  be  used  to  better  determine 
accelerometer  errors  and  improve  state  estimation.  These 
possibilities  were  not  explored  in  this  analysis  due  to  time 
limitations  and  a  feeling  that  confirmation  of  peak  vehicle 
speed  was  the  critical  area  of  concern. 


141 


Appendix  A.  Plots  from  Second  Iteration 
of  Extended  Kalman  Filter- 

This  section  presents  the  results  of  the  second  itera¬ 
tion  of  the  extended  Kalman  filter  developed  in  Chapter  III 
and  IV.  State  estimates  and  error  variances  at  the  initial 
time  from  the  first  iteration  of  the  Meditch  (5)  smoothing 
algorithm  are  used  to  update  the  initial  conditions  of  the 
extended  Kalman  filter  for  this  run.  In  addition,  a  slight 
correction  to  the  assumed  test  track  heading  provides  closer 
filter  correlation  to  range  and  azimuth  measurements,  allow¬ 
ing  a  reduction  in  the  estimated  measurement  noise  for  radar 
range. 

Figure  A.l  shows  a  plot  of  estimated  and  actual  radar 
(  range  obtained  with  the  corrected  test  track  heading.  This 

figure  indicates  that  the  divergence  of  actual  and  estimated 
range  between  20  to  24  seconds  has  been  removed. 

The  extended  Kalman  filter  shows  improved  convergence 
of  the  standard  deviations  of  the  state  estimates  due  to  im¬ 
proved  initial  conditions  from  the  smoother.  This  can  be 
seen  by  comparing  to  part  (b)  of  Figs.  A. 2  through  A. 8  to 
part  (b)  of  Figs.  4.17  to  4.23.  Figures  A. 9  (a)  and  A. 9  (b) 
are  plots  of  the  range  and  azimuth  measurement  residuals 
bracketed  by  the  residual  standard  deviations.  From  Fig. 

A. 9  (a)  it  is  apparent  when  the  residual  monitoring  routine 
bypasses  range  measurements  in  excess  of  three  times  the  re¬ 
sidual  standard  deviation.  Figure  A. 10  (b)  is  an  expanded 
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version  of  Fig.  A. 10  (a)  which  is  the  scalar  speed  estimate 
(magnitude  of  the  velocity  vector)  converted  to  Mach  number 
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Pig.  A.l(a)f(b)  Residual  range  analysis  after  correcting  track 
heading  to  179.209  degrees  true. 
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Pig.  A. 6(a), (b)  Extended  Kalman  filter  state  estimate  and 
standard  deviation  after  two  iterations 
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Pig.  A.9(a)f(b)  Range  and  Azimuth  residuals  and  standard 
deviations  after  two  iterations. 
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Fig.  A. 10(a), (b)  Velocity  converted  to  Mach  number  after 

two  iterations  of  the  Extended  Kalman  filter 


Appendix  B.  Optimal  Smoother  Computer  Program 


This  appendix  includes  a  listing  of  the  computer  pro¬ 
gram  used  to  incorporate  the  Meditch  (5)  backward-recursive 
optimal  linear  smoother  algorithm  described  in  detail  in 
Chapter  II.  The  data  required  for  this  program  includes  the 
current  sample  time,  t^,  state  estimate  vector  before  update, 
covariance  matrix  before  update  (stored  in  upper  triangular 
form) ,  state  estimate  vector  after  update,  and  covariance 
matrix  after  update  (stored  in  upper  triangular  form)  from 
tQ  to  tf. 

For  this  analysis  the  state  estimates  and  covariance 
before  and  after  each  measurement  update  are  stored  as  a  re¬ 
sult  of  one  run  of  the  extended  Kalman  filter  implemented  in 
SOFE.  A  short  data  reformating  program  was  used  to  put  the 
data  at  the  final  time  first  and  the  remaining  data  records 
in  backward-recursive  form  to  the  initial  time.  This  enables 
the  smoother  program  to  read  forward  through  the  data  but 

actually  compute  quantities  "backward"  in  time.  The  Meditch 

a  +  A 

algorithm  requires  x(t^  )  and  x(t^+1  )  for  one  calculation 
at  each  sample  time.  The  data  was  put  in  the  following  order 
to  allow  for  a  step-by-step  "read"  of  the  required  quantities 
for  each  time  of  calculation.  Each  data  record  contains  the 
following  information  in  the  order  shown: 

1.  time,  t^ 

_  A  + 

2.  state  vector  at  time  t^  af ter  update,  x(t^  ) 
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3.  covariance  matrix  at  time  t^  after  update,  P(t^+) 

A  _ 

4.  state  vector  at  time  t^  before  update,  x(t^  ) 

5.  covariance  matrix  at  time  t^  before  update,  P(t^-) 
The  order  of  the  system  used  in  the  rocket  car  analysis  al¬ 
lowed  for  storage  of  formated  records.  A  higher  order  sys¬ 
tem  will  probably  require  more  efficient  data  storage  and 
improved  formating  of  smoother  printed  output. 
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SUBROUTINE  MAKEP(N, PIN, POUT) 


COMPUTES  THE  “SMOOTHED"  COVARIANCE  MATRIX  AS 
IT  PROPAGATES  BACKWARD  IN  TIME. 

REQUIRES  MATRIX  ADDITION  AND  MULTIPLICATION 
ROUTINES  TO  BE  USED. 
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